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ABSTRACT 


An  analytical  study  is  conducted  to  determine  the  stability  and  free 
vibration  characteristics  of  laminated  anisotropic  plates  using  the  Levy 
approach.  Included  in  the  plate  model  are  the  effects  of  shear  deformation 
and  rotary  inertia.  Six  different  boundary  conditions  in  the  y  direction  are 
analyzed  in  oor junction  with  simply-supported  boundaries  in  the  x  direction. 
The  y  directed  boundaries  considered  are  simple-simple ,  clamped-clamped , 
sirnple-clamped,  simple-free,  clamped-free  and,  free-free. 

Solutions  are  presented  for  the  buckling  loads  and  natural  frequencies 
of  rectangular,  graphite-epoxy  symmetric  plates.  The  results  indicate  the 
importance  of  including  shear  effects  and  rotary  inertia  in  a  plate's 
mathematical  model.  The  overall  importance  of  these  equation  parameters  is 
definitely  a  function  of  the  boundary  condition  and  a  general  statement  cannot 
be  made.  In  addition,  the  effectiveness  of  the  Levy  technique,  in  studying 
laminated  problems,  becomes  apparent  in  handling  the  more  complicated 
boundaries  as  compared  to  the  Galerkin  or  Rayleigh-Ritz  techniques. 


M 


I.  Introduction 


Background 

Hie  use  of  ccrposite  material  in  ma rsy  varied  industries  has  increased 
tremendously  in  the  past  several  years.  This  is  largely  due  to  the  high 
strength  to  weight  ratios  of  ccnposites  as  well  as  their  ability  to  be 
tailored  to  meet  design  requirements  of  strength  and  stiffness.  A  special 
interest  is  the  use  of  ccrposite  materials  in  aircraft  structures. 

Coinciding  with  these  new  applications  is  the  need  to  better  understand 
the  physical  and  dynamic  responses  of  the  ccnposites  to  complex  in-plane 
stress  systems  [7],  Bast  research  has  clearly  indicated  the  need  for  a 
refinement  of  the  classical  plate  theory  (CFT)  in  order  to  better  predict 
ccrposite  plate  behaviour.  The  assunption  that  plane  sections  remain  plane 
after  deformation  (Kirchhoff  hypothesis)  results  in  a  mathematical  model  of 
plate  behaviour  which  is  overstiff  [7] .  The  need  to  include  through-the- 
thickness  shear  effects  was  first  recognized  by  Reissner  [14].  Mirdlin, 
shortly  thereafter,  added  the  thickness -dependent  effects  of  rotary  inertia, 
for  the  vibration  problem,  in  his  mathematical  model  for  the  flexual  motion  of 
isotropic  plates  [10].  Mindlin's  two  dimensional  theory  is  based  on  the 
premise  that  plate  displacement  is  a  result  of  two  rotations  due  to  bending 
and  two  rotations  due  to  shear  deformation.  Furthermore,  no  warping  of  the 
plane  sections  is  assumed.  This  inconsistency  is  somewhat  corrected  by  the 
introduction  of  a  correction  factor  [2]. 


Yang,  Norris,  and  Stavsky  [1]  used  Hindi  in's  model  to  develop  the 
frequency  equations  for  the  propagation  of  harmonic  waves  in  an  infinite,  two- 
layer  isotropic  plate.  Their  theory  (referred  to  as  YNS  theory)  was  applied 
by  Whitney  and  Pagano  [12]  to  laminated  plates  consisting  of  an  arbitrary 
number  of  bended  anisotropic  layers,  each  having  one  plane  of  material 
symmetry  parallel  to  the  central  plane  of  the  plate.  Their  study  [12] 
centered  on  solving  for  the  vibrations  of  antisymmetric  angle-ply  plate 
strips.  Work  by  Brunei le  [1]  on  transversely  isotropic  Mindlin  plates 
considered  the  stability  of  rectangular  plates  sinply-supported  on  a  pair  of 
opposite  edges  and  carrying  uniaxial  membrane  ocxipression  [24] .  This  work  is 
one  of  the  earliest  to  consider  the  feasibility  of  applying  Levy's  technique 
[13]  to  composite  plates  and  thereby  obtain  a  solution  for  displacement  that 
is  not  modelled  as  a  double  series  expansion.  The  accuracy  of  the  expression 
for  displacement  is  not  dependent  on  the  number  of  terms  retained  from  the 
expansion  in  this  case.  This  type  of  solution  is  ccmonly  referred  to  as  a 
closed-form  solution. 

Refinement  of  finite  element  analysis  with  the  inclusion  of  transverse 
shear  effects  began  with  the  work  of  Pryor  and  Barker  [25] .  They  studied  the 
cylindrical  bending  of  symmetric  cross-ply  laminates  using  a  model  based  on 
Reissner  plate  theory.  YNS  theory  was  modelled  by  Reddy  [18]  in  his  study  of 
the  free  vibrations  of  antisymmetric  angle-ply  plates.  His  later  work  [18] 
considered  orthotropic  laminates  of  bimodulus  materials.  A  different  approach 
was  used  by  Sathyamoorthy  and  Chia  [20]  in  their  study  of  anisotropic  skew 
plates.  They  applied  Von  Karman's  non-linear  plate  equations  to  calculate  the 
large  amplitude  vibrations  of  the  skew  plates. 


More  recent  work  by  Dawe  and  Craig  [6,7,16]  has  investigated  the  effects 


of  shear  deformation  in  a  number  of  plate  vibration  and  stability  problems. 

In  each  case,  the  Rayleigh-Ritz  or  finite-strip  methods  were  used  to  generate 
numerical  results.  Bowlus  [2]  analyzed  the  vibration  characteristics  of 
anisotropic  laminated  plates,  with  shear  deformation  and  rotary  inertia,  using 
the  Galerkin  method.  Reddy's  latest  work  [26]  applies  the  levy  technique  to 
the  bending  problem  of  symmetrical ly  laminated  rectangular  plates.  His  model 
includes  shear  deformation  and  considers  two  different  plates,  a  single 
layered  orthotropic  plate,  and  a  three  layered  cross  ply  composite  plate.  His 
numerical  results  are  generated  by  a  solution  based  on  the  state-space  concept 
developed  by  Franklin  in  1968. 

Thus,  the  literature  does  not  indicate  any  closed  form  solutions  for  the 
in-layered  synmetric  laminates  where  both  shear  deformation  and  rotary  inertia 
effects  have  been  considered  with  application  to  plate  buckling  and  plate 
vibration  respectively. 

Objectives 

There  are  three  main  objectives  to  this  thesis.  First,  the 
effectiveness  of  the  Levy  technique  in  calculating  natural  frequencies  and 
buckling  load  for  anisotropic  laminated  plates  is  determined.  The  plate's 
mathematical  model  includes  shear  deformation  through  the  thickness  and  rotary 
inertia.  Second,  the  technique  is  used  to  evaluate  "baseline"  solutions  for 
same  of  the  boundary  conditions  which  are  extremely  difficult  to  analyze  using 
the  approximate  methods  based  on  energy  principles.  Finally,  comparison  with 
classical  solutions,  when  available,  are  used  to  determine  the  importance  of 
shear  effects  and  rotary  inertia  in  the  plate  model. 


A  direct  approach  is  followed  to  successfully  realize  the  stated 
objectives.  The  notion  of  the  laminated  plate  is  modelled  using  YNS  theory 
[11].  Solution  of  the  derived  coupled  partial  differential  equations  of 
motion  is  obtained  by  application  of  the  Levy  technique.  The  plate  studied  in 
this  thesis,  a  specially  -  orthotropic  laminate,  does  not  contain  bending- 
extensional  or  bending-twisting  coupling  terms.  From  the  equations  of  motion, 
displacement  functions  cure  evaluated  and  used  to  solve  the  boundary  value 
problem  (BVP)  defined  by  the  application  of  specific  boundary  conditions. 
Simplification  of  the  BVP  leads  to  a  transcendental  equation  for  either  the 
natural  frequency  or  the  buckling  load,  depending  on  the  problem  under 
consideration.  A  computer  program  is  written,  in  Fortran  77,  to  solve  the 
transcendental  equation  using  the  incremental  search  method  [5] .  Limits  to 
the  effectiveness  of  the  Levy  technique  are  then  evaluated  in  terms  of 
allowable  laminates  and  plate  gecnetries/size . 

The  boundary  conditions  always  have  a  set  of  simple  supports  in  one 
direction  with  the  other  direction  consisting  of  either  simple-simple, 
clamped-clanped,  simple-clamped,  simple-free,  clamped-free,  or  free-free. 
Comparisons  of  results  to  classical  solutions  and  solutions  obtained  using  the 
various  energy  techniques  are  presented  whenever  available.  The  importance  of 
shear  deformation  and  rotary  inertia  is  determined  by  calculating  the  natural 
frequencies  and  buckling  loads  over  a  range  of  length  to  width  and  length  to 
thickness  ratios.  Initially,  rotary  inertia  is  neglected  in  order  to  better 
understand  the  relative  significance  of  considering  shear  deformation  effects. 
Rotary  inertia  is  then  re- introduced  into  the  vibration  plate  model  and 


calculations  are  repeated,  thereby  providing  an  indication  of  its  inportance 
in  the  accurate  modelling  of  thick  plates. 


II.  Theory  and  Modelling 

Anisotropic  Thick  Plate  Theory 

Classical  Laminated  Plate  Theory  (CUT)  incorporates  constitutive 
relationships  for  an  orthotropic  lamina  through  the  plate  thickness  resulting 
in  expressions  which  approximate  force  resultants  in  terms  of  displacement 
functions.  This  theory  provides  concepts  which  are  required  in  the  subsequent 
development  of  the  equations  of  motion.  As  a  starting  point,  the  basic 
constitutive  relationships  for  an  individual  lamina  are  developed  [2].  One 
should  refer  to  reference  [17]  for  a  more  conplete  and  detailed  derivation  of 
these  relations. 
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where  «1#  «2,  and  are  the  normal  strains,  «4,  «5,  and  eg  are  the  shearing 
strains,  Oy  a^,  and  cure  the  normal  stresses  and  a a^,  and  a&  are  the 
shearing  stresses.  The  terms  are  compliance  terms  and  may  be  written  in 
terms  of  the  lamina  engineering  contants  as: 
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in  the  ith  direction,  v 

.  is  Poisson's 

transverse  strain  in  the  jth  direction  when  loaded  in  the  ith  direction,  and 
Gjj  is  the  shear  modulus  in  the  i-j  plane. 

Equation  (1)  may  be  inverted  to  give  the  relationship  of  stresses  in 
terms  of  strain  in  the  form 

<*>-[Q']<<>  (3) 
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[Q'  ]  is  referred  to  as  the  reduced  stiffness  matrix  and  has  the  following 
form: 
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where 
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If  the  lamina  is  not  oriented  with  the  principal  x-y  axis  but  rather  is 
at  an  angle  e  (see  figure  2.1),  the  reduced  stiffness  matrix  must  be 
transformed.  The  matrix  applied  to  the  stiffness  terms  to  reflect  the  shift 
in  the  laminae  axes  is  defined  as: 


[T] 
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(6) 


where  1  =  cos  e  and  p  =  sin  9.  The  transformed  stiffness  matrix  for  the 
lamina  is  represented  as 


[Q'3  =  [THQ'ltT]1 
and  {a)  = 
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In  order  to  simplify  the  stress-strain  relationships  given  by  Eq(8) ,  a  state 

of  plane  stress  is  assumed  for  the  laminae.  That  is,  the  individual  laminae 

are  considered  to  be  thin  enough  that  the  average  value  of  a  across  the 

z 

thickness  is  negligible  [21].  Thus,  from  Eq.  (8) 

°z'0‘  °13'  ‘x  +  ®23'  *y  +  °33'  ‘l  +  ®36‘  \cy 
°r  .2  -  (Q13VQ33').x  ♦  <0J3*/6J3').y+  <^«V®i3'>>y  (10) 
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If  this  expression  for  «  is  substituted  back  into  Eq  (8) , 
the  resulting  transformed  stiffness  matrix  is  defined  as  [Q]  and 


has  the  form 
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With  the  removed  of  a  ,  the  stiffness  matris  becomes  a  5x5  matrix.  The 

z 

stress-strain  relationships  defined  in  Eq.  8  may  be  re-written  in  terms  of 
normal  and  shear  values  as: 
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These  relationships  are  for  a  single  lamina  and  may  be  used  as  a 

starting  point  in  the  derivation  of  expressions  for  the  forces  and  moments 

existing  in  a  laminate  of  N  perfectly  bonded  lamina.  Before  doing  so,  a 

further  simplification,  deeding  with  e_  is  discussed  as  well  as  the  convention 

z 


followed  in  defining  the  rotations 


and  ^y. 


The  strain  in  the  z -direction  is  assumed  to  be  small  enough  as  to  be 
negligible.  This  inconsistency,  widely  accepted  in  plate  theory,  implies  that 
no  stretching  occurs  in  the  direction  perpendicular  to  the  plate  midplane. 

For  a  laminate,  this  assunption  leads  to  a  discontinuity  in  e  at  the  lamina 
upper  and  lower  boundaries,  but  it  is  very  small.  This  assunption  allows  the 
modelling  of  the  displacement  field  for  the  laminate  with  YNS  theory.  Thus, 
the  plate  displacements  can  be  assumed  to  have  the  forms 

u  =  u°(x,y,t)+z*x(x,y,t) 
v  «  v°(x,y,t)+z0y(x,y,t)  (16) 


v  -  w(x,y,t) 


where  u,v,w  are  the  x,y,  and  z  coordinate  displacements  respectively,  u 
and  v°  are  the  pre-stressed  displacements  of  the  laminate  midplane,  and  and 
V>y  are  rotations  of  lines  perpendicular  to  the  midplane  due  to  bending  [2]. 

It  should  be  noted  that  the  inconsistency  mentioned  earlier  could  be  removed 
by  including  terms,  which  are  linear  and  quadratic  in  z,  to  the  expression  for 
w.  Reference  [9]  mates  it  clear  however,  that  for  most  plates,  the 
inconsistency  is  negligible. 

The  expressions  i/>  and  are  defined  as  rotations  about  the  y  and  x 
x  y 

axes  respectively  due  to  bending  moments.  With  the  axis  system  defined  in 
figure  2.1,  the  rotations  indicated  are  produced  by  positive  bending  moments. 


Figure  2.2.  Coordinate  System  of  Plate 


For  small  strains,  the  first  order  engineering  definitions  are: 


Substituting  the  displacement  field  given  by  Eg.  (16)  gives 
‘x  ”  u°'x  + 

‘y  '  v°'y  +  wy,y  <18) 

o  o 

7  suf  +  v  .  +  +  ztf 

xy  'y  'x  ^x,y  ^y,z 


vr-y 


+  Z  J  X. 


vhere  the  strains  at  the  plate  middle  surface  are 


<6y|  “ 1  v°'y 

7°xyl  u°'y  +  v°', 


and  the  midplane  curvatures  due  to  bending  are 


*y  *y,y 

rx,y  +  *y,x 


Egs.  (20)  and  (21)  give  expressions  for  the  strains  in  terms  of 
displacements  and  these  may  now  be  substituted  into  Eg.  (14)  to 
give  the  stress  -  displacement  relationships  for  the  kth  layer  of 
the  laminate: 


Q11  °12  °16 
®12  ®22  ®26 
®16  ®26  ®66 


66  k 


o  o 
U  ,y+V  , 


+Z  <  V>. 


-W> 

x.v  v.x 


4»:mV 


Following  Jones  [17] ,  the  resultant  forces  and  moments  acting  on  the 
laminate  are  obtained  by  integrating  the  stresses  in  each  layer  through  the 
laminate  thickness.  Thus,  for  a  laminate  described  in  figure  2.3,  with  N 
lamina,  the  forces  and  moments  cure: 


* 


Figure  2.3.  Geometry  of  laminate  consisting  of  N  lamina 
h/2  N  \ 

(Nx'W  =  f  (Wrxy)dZ  =E  1  (ax'Vrxy)dZ 

-ty2  tel  zk_1 

and  (23) 

W2  v  \ 

<*x'W  “J*  (WTxy)ZdZ  =  S  f  (VVrx y)zdz 

-h/2  lo=l  zk_1 

These  integrations  may  be  simplified  somewhat  as  the  stiffnesses  are 
constant  within  the  laminae  and  the  stresses  are  not  functions  of  z.  Thus, 
Eg.  (23)  my  be  rewritten  as 


*11  *12  *16 
*12  *22  *26 
*16  *26  *66 


B11  B12  B16 
+  B12  B22  B26 
B16  B26  B66 


B11  ®12  B16 
B12  B22  B26 
B16  ®26  B66 


D11  °12  °16 
+  D12  °22  °26 
°16  D26  D66 


where 


Aij  *  £  (0ij>K  'VV'l 


V2  2  (Qij>K 


1/3  £  (Qij>K  (ZX3-Vl3> 


Required  new  are  the  expressions  for  the  shear  forces  on  the  plate  in 
terras  of  displacements.  In  classical  plate  theory,  shear  deformation  through 
the  thickness  is  neglected  according  to  Kirchhoff  hypothesis.  In  this  thesis 
hewever,  Mindlin  plate  theory  is  used  which  allows  rotations  of  lines 
originaly  perpendicular  to  the  midplane  but  having  no  warping.  The  exclusion 
of  warping  is  incorrect  as  shear  varies  parabolically  through  the  thickness  of 
the  plate.  The  error  induced  by  this  inconsistency  is  reduced  to  an 
acceptable  level  by  the  introduction  of  a  shear  correction  factor  k.  The 
magnitude  of  k  was  calculated  by  Whitney  [15]  and  the  resulting  value  of  5/6 
is  used  in  the  constitutive  relations  for  transverse  shear. 


Fran  the  definition  of  engineering  strain,  the  relationships  for  the 


interlaminar  shear  strains  sure 


7W  =  3v  +  3w  =  +  w, 

^  3Z  8y  y 


7 


xz 


=  au  +  aw  =  v>  +  w, 
az  ax 


(29) 


(30) 


Substituting  Egs.  (29)  and  (30)  into  Eq.  (15) ,  and  introducing  the  shear 
correction  factor  k,  the  shear  stress  can  be  expressed  as 


-  - 

Tyz 

>=  k 

°44  Q45 

w'y  + 

rxz 

k  4 

®45  ®55 

v'x  *  *x 

The  resultant  shear  forces  through  the  thickness,  and  Q^,  are 
evaluated  by  integrating  the  shear  stresses  in  each  layer  through  the  laminate 
thickness.  Thus, 


h/2 

N 

*k 

S  Txz 

dz  =  Z 

J’  {Txz}k 62 

(32) 

-h/2 

k=l 

Vi 

h/2 

N 

*k 

S  ryz 

dz  =  Z 

f  {ryz>k  ^ 

(33) 

-h/2 

k=l 

Vi 

As  was  dene  for  in-plane  forces  and  moments,  the  integrations 


may  be  simplified  by  noting  that  the  stiffnesses,  Q. .,  are 


constant  within  the  laminae  and  the  shear  forces  are  independent 
of  z.  Thus 


°y  =  k  A44  A45  w'y  + 

°x  A45  ^5  w'x  + 

where  is  defined  by  Eq.  (26) 


(34) 


The  equations  of  motion  for  the  laminated  plate  may  now  be  devrived  as 

all  terms  appearing  in  those  equations  are  now  defined  in  terms  of  the  three 

unknown  displacements,  w(x,y,t) ,  V (x,y,t) ,  and  ^  (x,y,t) .  A  Newtonian 

x  y 

approach  is  used  to  derive  the  three  governing  equations  of  motion.  Figure 
3.4  presents  the  nomenclature  for  the  moment  and  transverse  shear  resultants. 


Figure  2.4.  Moment  and  Shear  Resultants  on  Plate 


Considering  the  forces  in  the  z -direct ion, 


s  Fz  •  'h  w'tt  ■  V^.x  ^-VW.y 


(35) 


(36) 


| 


& 


-  Nx  W'xx  +  Hy  M'yy 

p  =  plate  mass  density 
h  *  plate  thickness 


Taking  manents  about  the  x-axis 

-  Wy.tt  -  Wy  <V  -  Hy-Sy  dy  -  fiy>y  (dy)2 
+  V  +  V  <*  -  ”xy 
Taking  manents  about  the  y-axis 

*V  =  ^x,tt  “  *x  +  *x,x  *  -  \  -  v  +  v  +  Vy  * 


(37) 


“  °X  ^  "  °x,x 


W  ,  - 

where  I  =  J  pz  dz  =  ph  /12 
-h/2 


(38) 


(39) 


These  equations  may  be  simplified  if  higher  order  terms  are  neglected 


and  dx=dy=l  is  assumed.  Substituting  in  for  q,  Eqs. 

rewritten  in  the  form 

(35),  (37),  and  (38)  are 

°x,x  +  °y,y  +  NxW,xx  +  Ny  w'yy  =  ^'tt 

(40) 

*x,x  +  **xy,y  "  °x  =  T^y,tt 

(41) 

^xy^x  +  *y,y  -Qy  *  T*x,tt 

(42) 
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These  three  equations  can  new  be  expressed  in  terms  of  the  displacements 


by  substituting  Eqs.  (25)  and  (34)  into  the  equations  of  motion.  However, 
before  doing  so,  two  restrictions  must  be  applied  in  order  for  this  problem  to 
be  solved  using  the  Levy  technique.  This  thesis  considers  only  symmetric 
laminates  and  further,  only  those  which  are  special ly-orthotrcpic .  These 
restrictions  remove  all  coupling  stiffnesses  (B.^)  and  bending-twisting 
coupling  stiffness,  D16  and  D26,  reducing  Eqs.  (25)  and  (34)  to 


(  \  r  \ 


Dn 

D12  0 

*x,x 

N 

i  = 

D12 

°22  0 

\ 

*y,y 

W 

0 

0  V 

,  **y  +  Vx 

(43) 


1 

°y 

=  k 

A44  0 

W,  + 

'y  *y 

°x 

i 

-  0 

w'x  +  *x 

,  i 

(44) 


Assuming  the  time  dependence  of  the  displacements  to  be  harmonic  then 
e1^  allows  the  separation  of  the  time  variable  out  of  the  equations  of 
motion.  Thus,  the  equations  of  motion  can  be  written  in  a  form  containing 
only  displacements: 


D11  ^x,xx  +  D12  ^y,xy  +  D66  ^x,yy  +  D66  ^y,xy  ”  *^55  ^w'x  +  ^ 


+  coIV>x  =  0 


(46) 


°12  ^X/Xy  +  D22  ^y,yy  +  D66  ^x,xy  +  °66  ^y,xx  ”  ^44  ^w'y  +  'V 


+  U  Uy  =  0 


(47) 


where  u>  is  the  plate  natural  frequency 

3 

p  is  the  plate  mass  density  lb/ in 


These  are  the  three  coupled  partial  differential  equations  of  motion 
which  are  solved  for  the  plate  displacements  with  the  use  of  the  Levy 
technique. 
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LEVY  TECHNIQUE 


The  partial  differential  equations  describing  the  plate  displacement 
have  been  formulated  and  must  now  be  solved.  The  technique  employed  in  this 
thesis  to  solve  Eqs  (45) ,  (46)  and  (47)  is  the  Levy  technique.  This  technique 
is  unlike  sane  of  the  approximate  energy  techniques  (Galerkin,  Rayleigh-Ritz , 
Navier)  used  previously  to  study  the  problem  of  ocrposite  plate  stability  and 
vibration.  The  Levy  technique  leads  to  a  more  exact  solution  as  it  removes 
the  errors  associated  with  series  representations  of  the  variable.  Each 
displacement  term  in  the  equations  of  motion  becomes  a  single  unknown  value, 
as  opposed  to  a  series  of  unknown  variables.  The  mechanics  of  the  Levy 
technique  are  presented  in  the  next  section,  as  applied  to  the  specially 
orthotropic  plate,  in  order  to  clarify  some  of  the  general  characteristics  of 
the  method  presented  here. 


One  of  the  major  restrictions  of  this  method  is  the  requirement  that 
boundaries  on  two  opposite  sides  of  the  plate  be  maintained  as  sinple 
throughout  the  analysis  so  that  the  description  of  the  plate  displacements  may 
have  the  forms 


w(x,y,t)  =  [w(x,y)]e^ijt  =  z  W(y)  sin  n»x  (48) 

n=l  a 


too  ■ 

^x(x,y,t)  =  Wx(x,y)  =  Z  *x(y)  cos  mx  elwt  (49) 

n=l  a 


,  .  «  .  . 

*  (x,y,t)  =  [V>  (x,y)  =  Z  *  (y)  sin  mx  e1^  (50) 


% 

•y*! 


where  u>  =  natural  frequency 

These  forms  for  the  displacement  satisfy  the  mathematical  requirements 
for  edges  which  are  sinply-supported  at  x=0  and  x=a  given  by  Eq.  (51)  and  are 
used  to  simplify  the  equations  of  motion.  Note  that  the  time  variable  may  be 
factored  out  of  the  expressions  for  displacement.  This  leaves  functions  in 
ally  x  and  y  to  be  substituted  into  the  equations  of  motion. 

Simply  supported  edge  at  x=x  w(x,y)=  V>x(x,y)=  tfy'  (x,y)=0  (51) 


where  the  prime  denotes  differentiation  w.r.t.y. 


laminated  Plate 


Midi  has  been  said  on  the  benefits  of  using  the  Levy  technique  to  solve 
the  stability  and  vibration  problems  for  a  ocnposite  plate.  At  this  point, 
the  actual  simplifications  of  the  equations  of  motion  are  presented.  The 
resulting  equations  are  then  re-written  in  a  form  that  is  useful  in  solving 
either  the  stability  or  the  vibration  problem. 

In  order  to  vise  the  displacements,  as  described  in  Eqs.  (48) ,  (49) ,  and 
(50)  in  the  equations  of  motion  (45-47) ,  the  following  derivatives  are 
required  (note:  all  2  are  over  n=l  to  «,  no  time  dependence  shown) . 


w,  =  £  n*  W  cos  n*x 
a  a 


w,  =  2  W  sin  nrx 
y  a 


w'xx  =  2 1 - 


7J 


sin  mx 
a 


w,  =  2  W"  sin  n»x 
™  a 


K  v  =  E  2*  006  I** 

y,x  a  y  a 


*y.y  "  2  V  sl" 


=  2/- 


n2n\t  sin  mx 

~)  a 


XV  =  2  ^  V  0065  22- 
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A 


•!«1 


*  *  2  -r*  *v  Bin  mx 

x'x  a  x  a 


xt  =  2  *  '  cos  mx 
x,y  x  a 


*x,»  =  "/-J 


-n2*^\*x  006  mx 

~J  a 


* YW  =  2  *v'  sin  mx 

*/xy  a  x  a 


Vyy  -  s  V  0=6  JSS 
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where  all  primed  terms  denote  differentiation  w.r.t.y. 
Thus,  from  Eq.  (45) , 


r^l|T  ^  *x(Y)  kA55W(y)  +  **44  V  (y)  +  ^44  *"{Y) 

9 


-nV  Nx  W(y)  +  Ny  W"(y)  +  ,hA(y)J  si 


sin  mx  =  0 


from  Eq.  (46) 


z  H 

IFll 


i'JLL  D11  *x(y)  +  *  °12  V(y)  +  °66  V(y)  +  *  °66  *y'(y) 
la  2 


1 


1 


f  * 


If 


$ 

*,i 

$ 


I 

iQ 

I 

i 


“k  A55*x(Y)  k  ^55  W<y)  *■>  I#x(y) 

9 

from  Eq.  (47) 


1  006  mx  =  0 
a 


E  r] 

mil 


,("!*  °12  *x'(y)  +  D22*y"(y)  "  ^  D66*x'(y)  ~  —  D66*y(Y) 
L|  a  a2 

-  kA44*y(y)  -  J^44w'  (y)  ♦  w2I*y(y)  j  sin  M  =  0  (! 


Thus,  for  the  equality  to  hold,  the  following  equations  wist  be  satisfied  for 
each  n 


From  Eq.  (53) 


IDS.  WSs  *XM  -  n  *  **55  W(y)  +  *A44  *y'(y)  +  kA44W"(y) 


-n2*2  Hx  W(y)  +  Ny  W"(y)  +  ph^2  W(y)  *  0 


From  Eq.  (54) 


-n2*2  Dn  *x(y)  +  ^  °12  *y'(y)  +  °66*x"(y)  +  ^  °66  V(y) 

a2 


-**55  *x(y)  -  m  YK55  W(y)  +  </l*x(y)  *  0 
d 


From  Eq.  (55) 

-o*  d12  v<y>  +  D22V(y)  "  -  D**V<y)  - n2*2  D**vy> 


66  X 


66  y' 


-kA44*y(y)  -  *A44w'  (y)  +  «  »y(y)  ~  0  (58) 

Thus,  the  application  of  the  Levy  technique  reduces  partial  differential 
equations  (in  x  and  y)  to  ordinary  differential  equations  in  y.  Eqs.  (56), 

(57) ,  and  (58)  nay  now  be  used  to  solve  the  stability  or  vibration  problem  for 


the  specially^rthotropic  laminate. 


n 


u’ 

% 


i'!v 


i 

>$ 


I.] 

$ 

% 


1 


•Si 


Stability  Problem 


In  the  solution  of  the  stability  problem,  basically  a  static  problem, 
inertia  terms  will  not  be  retained  in  the  equations  of  motion.  In  addition, 
this  thesis  will  only  consider  the  case  of  a  unaxial  compression  force.  Thus, 

Ny  =  0. 


To  solve  the  equations  of  motion,  the  classical  approach  of  assuming  a 
very  general  form  for  the  displacements  is  followed.  Hence,  the  displacements 
are  taken  to  be 


W(y)  =  A*e0y 
*x(y)  =  B*e0y 
*y(y)  =  C*e0y 

*  *  * 

where  A  ,  B  ,  and  C  are  constants  to  be  evaluated. 


Further,  a  notation  is  adopted  which  allows  the  equations  to 
be  written  in  a  more  compact  and  manageable  form.  Using 


°12  +  °66  e 


nr/a  =  N 


Substituting  the  assumed  forms  for  the  displacements  and  appropriate 
derivatives  and  using  the  given  notation,  the  equations  of  motion  may  be 
simplified  and  re-written  in  matrix  form  as 


Note:  t  contains  the  function  e 

r  is  the  buckling  force  quantity 
end  n  contains  the  number  of  modes 

Ihe  norrtrival  solution  to  this  set  of  coupled  ordinary  differential 
equations  is  found  by  considering 

det  [Ay]  =  0  (62) 

At  this  point,  MACSYMA  [3]  is  used  to  calculate  the  factored  expression 
for  the  determinant.  The  Macsyma  package  is  available  on  the  AFITNET  and  a 
description  of  the  ocmnands  used  to  obtain  the  final  form  of  the  expression 
described  by  Eq.  (62)  is  presented  in  Appendix  A. 

The  resulting  equation  for  the  determinant  of  [Ay  ]  is  a  sixth  order 
equation  in  the  unknown  0.  It  can  be  written  as 


(64) 


where  the  coefficients  are  given  in  Appendix  A. 

This  can  be  reduced  to  a  cubic  with  the  substitution  of 


Thus, 

A3  f3  +  A2  f2  +  A1  f  +  Aq  =  0  (65) 

Eq.  (65)  can  be  solved  directly  for  three  roots  following  the 
trigonometric  technique  presented  by  Dickson  [4] .  The  roots  may  be  all  real 
numbers  or  a  single  real  number  and  two  complex  conjugates.  A  development  of 
the  problem  involving  three  real  roots  is  presented  first,  followed  by  a 
discussion  of  the  case  of  two  complex  conjigate  and  one  real  root. 

Under  the  premise  of  having  three  real  roots  to  Eq.  (65) ,  the  assumption 
is  made  that  two  of  the  real  roots  are  positive  and  the  third,  a  negative  real 
number.  Any  changes  to  this  assumption  are  discussed  in  the  section 
formulating  transcendental  equations  for  the  boundary  conditions. 

With  these  assumptions,  the  following  nomenclature  is  utilized 

(66) 

The  positive  roots  for  e  are  thus  assumed  to  be  a  and  7.  The 
displacements  described  in  Eq.  (59)  my  now  be  expressed  in  an  expanded  form 


„  Q  2  2 
ri  =  0l,2  =Q 

c2  =  e3,42=  (i^)2  =  -^2 


W(y)  »  A^e  J  f  +  Af  +  A^'-1  +  Age 

^(y)  =  Bxeay  +  B2e"13y  +  B3e^y  +  B4e_i^y  +  B5e7y  +  B6e~lY  (67) 

'l'y(y)  =  CxeaY  +  C2B~°y  +  C2e^Y  +  C4e-i^y  +  C5e7y  + 


As  presented,  Eq.  (67)  requires  the  evaluation  of  18  separate  constants 
for  each  set  of  boundaries  in  y.  Unfortunately,  the  bc's  define  only  six 
relationships  thereby  making  it  inpossible  to  evaluate  any  problem  using 
displacements  represented  in  the  form  of  Eq.  (67) .  However,  if  relationships 

•k  k  k  k 

between  A  /B  and  A  /C  are  established,  the  number  of  constants  to  be 
evaluated  drops  to  six  and  the  problem  becomes  well-posed.  Eq.  (61)  is  used, 

k  k  k 

for  this  problem  as  it  presents  the  relationships  between  A  ,  B  ,  and  C  in 
the  three  differential  equations  of  motion.  Since  these  three  equations  are 

k  k  k 

coupled,  the  simultaneous  solution  of  any  two  for  B  and  C  in  terms  of  A 
will  automatically  satisfy  the  third. 


Solving  the  second  equation  from  Eq.  (61)  for  C*  in  terms  of  A*,  B*: 


zA*  -  (fe2  -  z-N2 c 


Substituting  this  value  for  C  into  the  third  Eq.  from  (61)  and 

* 

simplifying  for  B 


[ (h92  -  z-N2 f)Nz  -  Neze2] 


[(hB-z-JTf)  (fe-z-Ng)  +  JTe2©2] 


(70a) 


and  thus  Nz-(fe2-z-N2g)A (0) 
Nie 


or  A(e)A*  =  C*  (70b) 

It  must  be  realized  immediately  that  A  (e)  and  A  (6)  are  row 
vectors  as  the  value  used  for  0  can  have  one  of  three  values  a, 

0  or  7.  The  choice  of  which  0  to  use  in  Eq.  (69b)  and  (70b) 
depends  on  the  argument  of  the  term  multiplied  by  A  (6)  or  A  (6)  in 
the  expressions  for  displacements.  Thus,  A (0)  is  equal  to  (a (a), 
A  (/9),  A(7))t.  The  discussion  is  not  yet  complete  as  a  potential 
problem  exists  when  calculating  A  09) . 

Since 

ip  =  0 

evaluation  of  A  09)  using  Eq.  (70a)  gives  an  expression  containing 
an  imaginary  number  in  the  denominator.  An  alternate  approach 
may  be  used  to  solve  for  the  expressions  of  A  (0)  and  a  09)  that  do 
not  contain  the  imaginary  number.  For  this  case 

*  .  *  .  * 

W(y)=A  suv9y  *  (y)=B  sip9y  and  *  (y)=C  cos^y  (71) 

x  y 

and  substitute  these  values  into  the  equations  of  motion.  Any 

two  of  the  coupled  equations  may  be  solved  for  the  relationships 

*  *  * 

between  A  ,  B  and  C  as  was  done  earlier.  The  calculations  lead 
to  the  simplified  expressions  for  A(0)  and  A(£): 

(h^+^f+z)  Nz-N ez02 

_  =  A  (/9) 


[N2e2^2-(tv32+N2f+z)  (f^+N2 gtz) ) 


(72) 
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Nz+(f^2+N2g+z)  AO) 


=  AO) 


The  displacements  in  terms  of  only  six  unknown  constants, 


are  respectively 


W  (y )  =A1eay+A2e^ay+A3ei^9y+A4e'‘i^y+A5e7y+A6e~r  Y 


*x(y)=A(a)  (A1eay+A2e‘ay)+AO)  (A3e^y+A4e'i^y)+A(7)  (A5e7y+A6e“Yy) 


*y(Y )=*(«)  (A1eay+A2e"Qy)+AO)  (A3ei^y+A4e'i^y)+A(7)  (A^+Age"77) 


With  the  trigonometric  and  hyperbolic  identities  of  Eq.  (75) , 


sinhx  =  1/2  (eX  -  e“x)  sinx  =  1/2 i  (e1*  -  e_lx) 


coshx  =  1/2  (ex  +  e_x)  cosx  =  1/2  (e1*  +  e-ix) 


the  displacement  equations  nay  be  written  in  a  more  recognizable 


form.  For  example,  the  displacement  in  the  z  direction,  W(y)  is 


?y  +  A.e“y  +  A.e-°y  +  a  e^y  +  a  e^y  +  a  e^y 


W(y)  =  A.e  J  +  A.e  J  +  A_ 
2X  2 1  2* 


A_e  +  ***  *  +  a  € 
2T  2J  2J 


+  A  e-;^y  +  A  e~^y  +  Aj.e7y  +  a  e7y  +  A  e_7y  +  A  e_7y 
«4  ^r5  .3  .6  .6 


or  alternately 


W(y)  =  (A1+A2)l(eay+e'iay)  +  (A1-A2)l(eQy-e"Tiy)  +  (A3+A4)l  (e^y4e_i^y) 


+  (A3-A4)i  (ei^y-e'i^y)+(A5+A6) l(e7y««"ry)+  (A5~A6)l(e7y-e“Ty) 


applying  the  identities  of  Eq.  (75) ,  W(y)  becomes 

W(y)  =  Asirhay+Booshay+<Six\8y+Dcx3S|i9y+Esirih7y+Fcx3Sh7y  (76) 

where  A=  A^-A^  B^+Aj,  OfA^A^i 
D=  A3+A4,  E^-A^  F^+Ag 

Similarly,  *x(y)  and  ty(y)  become 

*x(y)=A  (a)  (Asinhay+Booshay)+A(0)  (Csir\fly+Doos0y) 

(77) 

+A  (7 )  (Esinhyy+Fooshiy) 

*y(y)=fi  (a)  (Acoshxy+Bs inhay)  +A  (fi)  ( -Cocs^y+Dsir^y) 

+a(7) (Bcosh7y+Fsinh7y)  (78) 

It  is  important  to  note  that  W(y)  and  *  (y)  are  even  functions  and  *  (y) 

x  y 

odd.  This  is  dictated  by  the  number  of  derivatives  of  each  expression  defined 
in  the  equations  of  met  ion.  Eqs.  (76) ,  (77) ,  and  (78)  represent  forms  for  the 

displacements  which  may  be  evaluated  for  specific  boundary  conditions  in  y. 

.. 

Eqs.  (76) ,  (77) ,  and  (78)  thus  represent  forms  for  the  displacements  which  may 
be  evaluated  for  specific  boundary  conditions  in  y.  These  forms  have  been 
derived  under  the  premise  that  all  three  roots  of  Eq.  (65)  are  real  numbers. 

For  the  case  where  the  roots  of  Eq.  (65)  are  complex  conjugates  and  a 
single  real  number,  the  development  of  corresponding  expressions  for  the 
displacements  is  not  straight  forward.  The  problem  lies  in  deriving  the 
proportionality  functions  relating  A*/B*  and  A */C*.  Without  these  functions, 


the  number  of  unknown  constants  which  must  be  evaluated  at  each  boundary 
condition  is  eighteen  and  the  problem  is  not  solvable.  Thus,  the  author  has 
recognized  the  possibility  that  specific  roots  nay  be  missed  under  the 
assumption  that  all  roots  to  Eq.  (65)  are  real  numbers.  The  overall  results 
when  plotted  however,  will  indicate  trends  for  a  number  of  different  values 
for  a  parametric  function.  The  potential  of  answers  not  conforming  to  the 
associated  technique  is  present  but  the  values  will  be  smeared  out  in  the 
extrapolation  of  the  function.  Thus,  the  overall  trend  is  extended  to  any 
area  of  difficulty  and  yields  results  which  are  within  tolerance  of  any 
obtained  by  trying  to  reformulate  the  problem  and  then  solving  for  the 
appropriate  functional  forms  of  the  displacements  when  only  one  real  root 
exists.  The  physical  nature  of  the  buckling  or  vibration  phenomena  studied  in 
this  thesis  indicates  a  lack  of  irregularities  in  the  solution  and  thus 
substantiates  this  approach  to  dealing  with  complex  conjugate  roots  to  Eq. 

(65) .  The  transcendental  equations  for  each  of  the  boundary  conditions 
studied  cue  presented  following  the  development  of  corresponding  forms  of  the 
displacements  for  the  vibration  problem. 


Vibration  Problem 


For  the  solution  of  the  vibration  problem,  the  procedures  followed  will 
parallel  exactly  that  presented  earlier  in  the  stability  analysis.  Inertia 
terms,  translational  and  rotary,  will  be  retained  in  the  equations  of  motion. 
Hie  in-plane  load  Nx  will  of  course  be  excluded  from  the  derivation  since  the 
lowest  natural  frequency  will  be  affected  by  its  presence  in  the  equations  of 
motion. 


Hie  forms  chosen  for  the  displacements ,  given  by  Eq.  (59)  are  used  to 
solve  the  equations  of  motion.  Additional  notation  is  introduced,  to 
supplement  that  presented  in  Eq.  (60) ,  representing  the  inertia  terms 


J2!  =  x 
.  2 

phw  =  y 


(79) 


Using  previously  defined  notation  and  Eq.  (79) ,  the  equations  of  motion 
in  matrix  form  for  the  vibration  problem  are 


y+t^-N^z 

-NZ 

tz 

A 

A 

-Nz 

x+ft2-gN2-z 

eNt 

* 

B 

-tz 

-eNt 

x+ht2-fN2 -z 

M 

(80) 
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Considering  the  nontrivial  solution  to  this  set  of  equations,  the 
expression  for  det[A^_.  ]  is  presented  in  Appendix  B.  The  resulting  sixth  order 
equation  in  e  may  be  reduced  to  a  third  order  equation  by  the  substitution 
defined  in  Eq.  (64) 

The  cubic  in  $■  is 


V3  +  b/  +  Blf  ♦  B0  =  0 


!i*2i 

5| 


where  the  coefficients  Bi  are  given  in  Appendix  B. 
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Jfjj 


The  roots  of  Eq.  (81)  can  be  expressed  in  terns  of  a,  p ,  and 
7  and  the  expanded  expressions  for  the  displacements  are 
identical  to  those  described  in  Eq.  (67) .  Once  again,  a 
relationship  is  sought  between  the  constants  A*,  B*  and  C*  and 
the  expressions  corresponding  to  A(0)  and  A(0)  for  the  vibration 
problem  are  given  by 


[he^z-N2  f+x)Nz  -  Neze2] 

[  (h^-z-l^f+x)  (f^-z-N2 g+xJ-HJ^e^2] 


*  * 
A  -  B 


(82a) 


n(e)A*=B* 


(82b) 


ia 

w 

& 


and  Nz-(fe2-z^i2g+x)n(e) 


-  *  * 
n(6)A  =  c 


*  * 
A  =  C 


(83a) 


(83b) 


As  before,  the  alternate  forms  required  for  the  case  of  e=ip 


are  given  by 


(l\92+N2f+z-x)Nz  -Nezp2 

_  =0(0)  (84) 

[N2e2^2-(h02+N2f+z-x) (f^2+N2 g+z-x) ] 

and 

Nz+  (N2  g+f02+z-x)  n  (p ) 

_  =  0(0)  (85) 

Hep 


The  final  forms  for  the  displacement,  shown  in  Eqs.  (76) ,  (77) ,  and 
(78) ,  remain  basically  the  same  the  only  difference  being  the  multiplication 
factors.  Thus, 


W(y)  =  Asinhay+Booshay+CSinSy+Dcos^y+Esinhry+Fooshyy  (86) 

#x(y)  =  Q(q)  (Asinhtay+Booshay)  +o  (p)  (CSirv3y+Dcos/3y) 

+  0(7)  (Esinhry+Fcx5sh7y)  (87) 

Ofy(y)  =  0(a)  (Aoosttoy+Bsinhay)  +0  (p )  (-Ocos0y+Dsir$y) 

+  0(7)  (Ecoshyy+Fsirirry)  (88) 


These  equations  are  used  in  establishing  the  transcendental  equations 


for  the  natural  frequencies. 


As  the  equations  for  the  displacements  are  similar  for  the 


stability  and  vibration  problems,  transcendental  equations  for 
each  boundary  condition  in  y  are  developed  only  for  the  stability 
case.  Corresponding  equations  for  the  vibration  analysis  follow 
directly  by  substituting  0(0)  for  A(0)  and  0(0)  forA(0).  Each 
boundary  condition  in  y  is  considered  independently  and  the 
simplified  equations  for  the  displacements,  in  terms  of  only 
three  constants,  are  presented  in  matrix  form.  The 
transcendental  equation  for  each  case  is  obtained  by  taking  the 
determinant  of  the  three  by  three  matrix. 

Simple-Simple  Boundaries 

For  the  case  of  simple-simple  boundary  conditions  in  y,  the  following 
physical  and  reesulting  mathematical  conditions  are  used  to  define  the 
constants  in  the  displacement  expressions. 

Physically  at  y=y  w(x,y)  =  V>x(x,y)  =  M^.(x,y)  =  0 

Mathematically  W(y)  =  *x(y)  =  #y(y)  =0  (89) 

Note:  y  represents  a  general  y  position  of  the  boundary 

If  Eq.  (89)  is  used  to  solve  for  the  constants  in  the  displacement 
equations,  the  following  results  are  obtained: 

A=B=0=E=F=0 

and  for  a  non- trivial  solution 
sin  fib  -  0 
or  fib  =  n* 


Thus,  for  the  s inple-s iraple  boundaries  in  y  there  two  positive,  one 
negative  real  roots  to  Eq.  (65) ,  the  relationship  that  is  programmed  and 


solved  for  is  straight  forward: 

det  =  p-  w/b 
and  n  =  1 


(90) 


If  three  positive  real  roots  exist,  the  expression  that  must 
be  solved  at  y=b  becomes 

W(b)  =  0  =  AsinhabKSirh^b+Esinhyb 
*x(b)  =  0  =  Av1sihh>bKV2sinh^bfEv3sinh7b  (91) 

*y(b)  =  0  =  Aav4sinh»bK^v5sinlT^b+E7V6sinh7b 

Thus,  the  sirv9y/oos0y  terms  in  the  displacements  are  replaced  by 
sinhfly/coshfly  if  the  third  root  is  positive  real. 


The  matrix  expression  for  Eq.  (91)  is 


[a^] 


{0} 


(92) 


where 

all 

=  sinhab 

a12 

a21 

=  v^sinhab 

a22 

an 

=  v.asinhab 

=  sinhab 

a13  =  sinhyb 

=  v2sinh9b 

a23  “  v3sinbyb 

=  v^sinhSb 

a  =  v  rsinhyb 

Boundaries 


The  following  physical  and  mathematical  conditions  define  a 
clanped  boundary  at  y=y 


Physically  w(x,y)  =  ^x(x,y)  =  i/>y(x,y)  =  0 

Mathematically  W(y)  =  *  (y)  =  *  (y)  =  0  (94) 

x  y 

Thus,  at  y=0,  the  displacements  become 

W(0)  =  0  =  BfDfF 

*x(0)  =  0  =  Bv1+Dv  +FV  (95) 

*  (0)  =  0  =  Av.-Cvc+Ev, 
y  4  d  o 


Solving  Eg.  (95)  in  terms  of  the  three  constants  yields 


V  —v 

D  «  B  1  3  =  Bv 

V3~V2 


v  —v 
2  1 

F  =  B  _ _ _  =  Bv, 

v  — v 
3  2 


8 


E  =  C  5  -  A  4 


(96) 


(97) 


(98) 


Using  the  relationships  described  by  Eq.  (96),  (97),  and  (98)  the  BC's 
equations  can  be  determined  as  a  function  of  A,  B,  and  C  at  y=b.  This  yields 
the  following  matrix  equation: 


where 


a11  =  sinhab  -  v4  sinhrb 


a12  =  ooshab  +  vyoos^b  +  VgOoshyb 


a13  =  sirv3b  +  sinhyb 


a21  =  v^inhab  -  v1  v4  sinh-fb 


a22  =  v1  ooshab  +  v2  v?  aos^b  +  v3vg  ooshrb 


a23  =  v2  sirv3b  +  v3  v5  sinhyb 


a31  ‘  V1  °0Sh“b  *  V4OOSh’b 


a32  =  v4  sinhab  +  v?  sirtfb  +  vg  vg  sinhrb 


ag3  =  -v5  ooefih  +  ooshyb 


where  a,  0  and  7  are  variables  and  are  functions  of  these 
variables 


Boundaries 


The  mathematical  conditions  describing  the  simple  and  clamped  boundaries 
are  given  by  Egs.  (89)  and  (94).  Considering  a  simple  boundary  at  y=0  leads 


\ y*v%y 


»-*  'r  ■'tri,' 


i  v4  >*  „ 
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i 
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W(0)  =  0  =  BfEM-F 

tx(°)  =  o  =  Bv1+Dv2+Fv3 

(100) 

*y' (0)  =  BaV4+qffv5+F7V6 

or 

09 
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II 

►*» 
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o 

(101) 

A  clamped  boundary  at  y=b  gives 

W(b)  =  0  =  Asir^bKSirvSbfEsinhrb 

*x(b)  =  0  =  Av1sinhab+Cv2sir>fibfEv3sinh7b  (102) 

♦y  (b)  =  0  =  Av4coshab-Cv5ccs^9b+Ev6cosh7b 

Eq.  (102)  in  matrix  form  becomes 

=  {o}  (103) 

where 

a1]L  =  sinhab  a^  =  sirvSb  a13  =  sinhyb 

a21  =  v^inhab  a22  »  v2sir>8b  a23  =  v^inh^b 

a31  =  v4ooshab  a32  =  -v^cos&b  a33  =  VgCosh^b 

Simple-Free  Boundaries 

A  siitple  boundary  at  y  =  0  leads  to  the  relationships  between  constants 
given  by  Eq.  (101) .  A  free  boundary,  at  y  =  b,  has  the  following  physical  and 
mathematical  conditions 


[a.  .] 
L  il 1 


< 

pi 

f: 

*• 

f 

* 


i 


Htysically 

Mathematically 


=  o 


JV  =  -ND12*x(b)+D22V(b) 


=  *x'(b)+»y(b)  =  0 
Qy  =  W(b)4*  (b)  =  o 
43 


=  0 


(104) 


Expanding  Eq.  (104)  and  placing  in  matrix  form  gives 


where 


'VIS)'  ,0) 

(D22aV4-ND12Vl)  sirftob 

(P22^V5“ND12V2,si^b 

(D227V6~ND12V3)  sinhT'b 
(qv1+Nv4)  coshab 

(/9v2^Nv5)  cos^b 

(7V3+Nvg)  cosh7b 

(a+V4)  coshab 

(^-v5)cos^b 

(7+Vg)  COSh7b 


(105) 


Boundaries 


A  clamped  boundary  at  y=0  gives  the  relationships  described 

by  Eqs.  (96) ,  (97) ,  and  (98) 

D  =  Bv_  (96) 

F  =  Bv '  (97) 

E  =  CVg/Vg  -  AV4/V6  (98) 

Substituting  these  relationships  into  the  expanded  form  of  Eq. 

(104)  and  rewriting  in  matrix  form  gives 


[aij]  jB  =  {0} 


(106) 


where 


all  -  <-™12VD22“V4>si"h°b  +  <v3  V4ND12 


D22tv4)  sinhrb 


in.  .  v ^  \  l*.  v rw  tJ-fe  vr* 


ai2  =  (-^l2VD22aV4)0C^bf(-ND12V2+D22^V5)V7C0^b 


+  (-ND12V3+D227V6)V8°°S^b 
a13  =  (■ND12V2+D2/V5)  (-NDi2v3v5+D22^v5)  si^b 

a21  =  (aV^"*"^V4  )  COShab-  (7V3V4+V4 )  OOSltyb 

a22  =  (av1+Nv4)sinhab-Ov2+Nv5)v7sir\9b+-(7V3+Nv6)v8sijnh7b 
a23  =  O?v2-Nv5)cos0bf  (7v3v5+JV5)oosh7b 

a31  =  (a+v4)cx>shQb-(7V4+v4)cxish7b 

a32  =  (Q+v4)sinhab+(-^+v5)v7sirv9tH-(7+v6)vgsinh7b 
a33  =  (^-v5)cos)9bf(7V5+v5)cosh7b 


Free-Free  Boundaries 


A  free  boundary,  as  defined  mathematical ly  by  Eq.  (104) 


gives  for  y=0 


My  =  ^JD12[v1B+-v2Dfv3F]+D22[v4QB+-v5^I>fv67F]  =  0 
=  v^A+v^C+v^E+N  [  v^A-VgO+VgE  ]  =  0 
Qy  =  v4A-v5C+VgE+Ai+Q9+E7  =  0 

from  Eq.  (107) ,  E  can  be  written  in  terms  of  A  and  C  as 
-(Nv,  +  v,)  (Nvc  -  v„) 


(v3  +  Nv6)  (v3  +  Nv6) 


't*  +  V  = 


(111) 


*2*1 

|0| 

IV 


With  this  relationship  and  Eq.  (109) ,  a  relationship  between 
A  and  C  is  developed 


(V4+a)  -  (Vg+7)  (NV4+Vx) 

A  + 

(£-v5)  +  (v6+-y)  (Nv5-v2) 

C  =  0 

(Nv6+v3) 

(Nv6+v3) 

(112) 

or 

V9A+V10C  =  0 

(113) 

thus 

C  = 

“V  -  V11A 
^0 

(114) 

and  it  follows  from  Eq.  (Ill)  that 

<v7  +  V8V. 


(V7  +  Vll*  A  =  E 


V12  A  “  E 


(115) 


from  Eq.  (107) 


<D22»VND12V1>  „  x  <D22^5-®12V2» 

15  ■+■ 


(ND12V3‘D227V6) 


(ND12V3‘D227V6) 


D  =  F 


(116) 


V13B  +  V14°  =  F 


(117) 


Introducing  the  relationships  between  the  constants  to  the  relationships 
for  a  free  boundary  at  y=b  and  expressing  in  matrix  form  leads  to 


[a.  .] 
L  ir 


{0} 


(118) 


where 

au  =  (-ND12v1+D22v4a)sinhQb+(-ND12v2v11+D22v5v11^)sin^b 
+  (-NDi2V3v12+D22V6V127)  Sll^b 

a!2  =  (-^12VD22V)C^bf^12Vl3+D22Vl37)0^b 

ai3  =  (“^12VD22V)^bf(-^12V3V14+D22Vl4^°^b 

a2i  =  (Vxa4-Nv4 )  ooshabf  (V2Vn'9_Nv5V1i)  (V3V127+NV6V12 )  oosl^b 

a22  =  (v1q+Nv4  )  sinhabt-  ( v3v13Tr+Nv6v13 )  sinhyb 

a23  =  ( -v^+NVg)  sinflbf  (  V3v147+Wv6v14  )  sinhrb 

a31  =  (v4^)oosfebf(-v5v11+^v11)°osh^bf(v6v1247v12)ooshrb 


a32  =  ( v44tl )  sinhab+  (v6v13+7V13  )  sinhrb 
a33  =  (v5-^)siJiSbf(v6v14+7v14)sinh7b 


III.  Discussion  and  Results 


This  chapter  presents  a  brief  introduction  of  the  two  corputer  programs 
written  to  solve  the  transcendental  equations,  generated  earlier,  for  either 
the  stability  or  vibration  problem  of  a  laminated  plate.  Discussions  follow 
on  the  physical  properties  of  the  plate  as  well  as  the  analysis  performed 
using  the  plate. 

Computer  Programs 

One  existing  computer  program  is  modified  and  one  program  is  written  to 
solve  the  respective  stability/vibration  problem  of  a  rectangular  plate  with 
the  six  different  boundary  conditions  already  discussed.  The  first  program  is 
a  modification  of  a  program  written  by  Bowlus  [2]  which  calculates  the 
nondimens ional  bending  and  extensional  stiffnesses  for  a  symmetric  laminate. 
The  second  program  formulates  the  boundary  value  problem  for  a  particular 
boundary  condition  in  y  in  the  form  [a^]{c}={0}  where  [a^  ]  is  the  3x3  matrix 
containing  the  eigenvalue  and  c  is  a  column  vector  of  displacement  constants. 
The  program  solves  the  transcendental  equation  given  by  det  [a  —  ]  for  the 
value  of  Nx  or  w.  Appendix  D  gives  a  more  detailed  description  of  the  second 


program. 


Analysis  Performed 


TVk>  main  areas  are  analyzed  in  this  thesis,  using  the  data  generated  by 
the  second  program.  First,  the  ability  of  the  Levy  technique  to  effectively 
solve  the  various  vibration  and  stability  problems  formulated  earlier. 

Second,  the  importance  of  shear  deformation  (and  for  the  vibration  problem 
rotary  inertia)  in  a  mathematical  model  of  a  laminated  plate. 

The  application  of  the  Levy  method  is  validated  for  the  vibration 
problem  by  comparing  results  to  those  calculated  by  Bowlus  [2] .  The  boundary 
conditions  used  in  this  comparison  will  be  simple-simple  in  the  y-direction. 
Due  to  the  similarities  in  the  formulation  of  the  vibration  and  stability 
problems,  such  a  comparison  will  also  validate  the  stability  portion  of  the 
program.  The  author  could  not  find  other  published  works  for  plates  with 
similar  gecmetries/material  properties  and/or  ply  lay-ups  and  different 
boundary  conditions  with  which  to  further  evaluate  the  second  program. 

General  trends  and  expectations  for  specific  boundary  conditions  are  available 
however,  and  are  used  in  the  later  discussions  of  the  results. 

The  impact  or  importance  of  shear  deformation  is  evaluated  by  altering 
the  length  to  thickness  ratio  of  the  plate.  For  the  vibration  problem,  rotary 
inertia  is  introduced  into  the  model  and  a  second  set  of  calculations 
performed.  For  three  of  the  boundaries,  notably  sinple-sinple,  simple-clamped 
and  simple-free,  a  square  plate  is  used  in  the  calculations.  For  the  other 
three  boundary  conditions  in  y,  the  formulation  of  the  transcendental 
equations  forces  the  use  of  rectangular  plates  with  aspect  ratios  of  two  or 
greater  due  to  computational  limitations.  The  analysis  follows  a  presentation 


of  the  material  properties  of  the  laminated  plates  used  in  this  analytical 
study. 


laminated  Ccrrposite  Plate  Properties 

The  plate  studied  in  this  thesis  is  constructed  of  a  graphite-epoxy 
(AS/3501)  material  and  has  the  following  material  properties 

E1  =  21.0E+06  psi 
E-  =  1.40E+06  psi 

(119) 

G12  =  0.6E+06  psi 
,12  =  0.3 
p  =0.055  lb/ in3 


The  plate  has  a  ply-layup  of  [0,90]^.  For  the  ccnparison  with  Bcwlus, 
m  equals  100  indicating  a  total  laminate  thickness  of  one  inch  if  each  lamina 
is  assumed  to  be  0.005  inch  thick.  A  plate  with  m  equal  to  200  is  used  for 
all  other  calculations.  Tables  3.1  and  3.2  contain  the  stiffness  values 
obtained  from  the  first  program  for  the  two  plates. 
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One  Inch  Thick 


Parameter 

Value 

A44 

540,000.3125 

^55 

540,000.3125 

D11 

1555164 . 375 

°12 

35211.2656 

°22 

322769.875 

°66 

50000.0 

Units  for  A^j  :  lb/ in 


Units  for  D^j  :  lb-in 


Table  3.1  Stiffness  Parameters  for  a  One  Inch 
Thick  Plate 


Graphite-Epoxy  [0/90] 


200  s 


Two  Inches  Thick 


Parameter 

Value 

A44 

1,080,000.625 

^55 

1,080,000.625 

D11 

12,441,315.0 

°12 

281,690.125 

°22 

2,582,159.0 

°66 

400,000.0 

Units 

for  Ay 

:  lb/ in 

Units 

for  Dy 

:  lb-in 

Table  3.2  Stiffness  Parameters  for  a  Two  inch 
Thick  Plate 


Characteristics  of  Levy 


In  most  of  the  previous  studies  of  laminated  plate  behavior,  proof  of 
convergence  of  the  solutions  must  be  presented.  With  the  Levy  technique,  the 
values  of  or  are  calculated  using  a  closed  form  algorithm.  The  only 
factor  that  nust  be  ensured  is  that  the  value  computed  is  indeed  the  lowest 
value  for  oru.  This  is  done  by  repeating  the  computations  for  increasing 
values  of  n,  the  bending  mode  number,  until  the  user  is  satisfied  that  the 
results  are  continuously  increasing.  The  lowest  eigenvalue,  regardless  of  n, 
is  the  sought  after  solution.  In  most  cases,  n  need  not  be  increased  past 
three  to  determine  the  trend  of  the  output. 

In  order  for  the  program,  which  is  used  to  evaluate  the  transcendental 
equation,  to  have  some  credibility,  results  must  be  compared  with  previously 
published  work.  As  was  stated  earlier,  apart  from  work  by  Bcwlus,  published 
work  on  the  buckling  or  vibration  of  [0/90]ms  composite  materials  was  not 
abundant.  Leissa  [27,28]  and  Brunelle  [1]  have  results  which  can  be  used  to 
validate  trends,  but  only  Bcwlus  has  numerical  values  which  can  be  used  in  a 
direct  comparison.  Thus,  from  reference  [2]  for  a  plate  simply  supported  on 
all  sides,  the  following  comparison  cam  be  made  for  the  first  natural 
frequency. 


s 

(a/h) 

Galerkin  Method 
with  SD  and  Ri 

Levy  Method 
with  SD  and  RI 

% 

Difference 

17.5 

3754.98 

3762.44 

0.20 

20.0 

2899.69 

2913.32 

0.47 

22.5 

2310.70 

2320.17 

0.41 

27.5 

1559.94 

1568.82 

0.57 

30.0 

1321.80 

1322.56 

0.06 

35.0 

975.16 

976.10 

0.10 

40.0 

752 . 80 

749.47 

0.44 

Table  3.3  Comparison  of  Galerkin  and  Levy  Techniques 
for  a  Plate  Simply-Supported  On  All  Four 
Edges.  Plate  Thickness  is  One  Inch. 

The  values  in  Table  3.3  shew  good  agreement  between  the  closed  form 
solution  and  one  obtained  using  the  Galerkin  technique  with  the  double  series 
containing  six  terms  each. 

To  further  validate  the  second  program,  a  comparison  is  presented  of  the 
results  computed  for  the  sinple-claiped  boundary  in  y  and  the  classical 
solution  obtained  for  the  specially  orthotropic  plate  using  an  equation  from 
Whitney  [8],  For  a  plate  with  no  shear  deformation  or  rotary  inertia,  the 
natural  frequency  cam  be  computed  from 


where 


R  =  a/b 
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=  (nH-.25)w 

a2  =  n2’r2a1(a1“l)  (121) 

=  I> 

a  =  80  in 
p  =  0.055  lb/ in3 

for  all  m  and  n.  (for  fundamental  frequency,  m=n=l) 


Using  the  properties  for  the  two  inch  thick  plate  given  in  Table  3.2  and 
computing  the  natural  frequency  for  a  plate  where  R=l.O,  the  following  is 
obtained 


u  =  769.574  Hz 


(122) 


This  value  represents  the  classical  plate  solution  and  can  be  ccnpared 
to  the  value  obtained  for  the  natural  frequency,  using  the  Levy  technique, 
given  in  Table  3.4 


s 

(a/h) 

Levy  Method 

SD  and  RI 

Levy  Method 
SD  no  RI 

10.0 

4891.466 

4917.118 

12.5 

3282.50 

3294. 

15.0 

2344.26 

2350.771 

17.5 

1753.1 

1756. 

20.0 

1358.306 

1360.588 

22.5 

1082. 

1083. 

25.0 

881.889 

882.872 

27.5 

732. 

732. 

All  frequencies  are  in  Hz 

Table  3.4  Natural  Frequency  of  Two  Inch  Thick  Plate 
With  and  Without  RI.  Sinple-Clanped 
Boundaries  in  Y. 

Thus,  as  the  plate  thickness  ratio  increases  and  the  effects  of  shear 
deformation  and  rotary  inertia  become  negligible,  the  results  approach  the 
classical  solution  of  769.574Hz.  Classical  solutions  for  other  boundary 
conditions  or  for  the  stability  problem  could  not  be  found  in  the  literature. 
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Analysis  of  Shear  Deformation  and  Rotary  Inertia 

The  second  nain  area  of  investigation  in  the  thesis  is  the  importance  of 
including  shear  deformation  in  the  plate  model  for  the  six  different  boundary 
conditions  proposed  for  the  plate  study.  In  the  case  of  the  vibration 
analysis,  the  benefits  of  futher  model  refinement,  due  to  the  addition  of 
rotary  inertia,  is  also  investigated.  The  analysis  is  accomplished  by 
performing  calculation  for  various  plate  length  to  thickness  ratios  and 
comparing  any  trends  to  behavior  predicted  using  classical  plate  theory. 

The  following  discussion  is  separated  into  two  sections;  the  first 
considers  the  buckling  problem,  the  second,  the  free  vibration  problem. 

Buckling  Problem 

For  each  boundary  condition,  increasing  length  to  thickness  ratios  and 
length  to  width  ratios  cure  programmed  and  solved.  Computational  limitations 
play  a  limiting  role  as  feu:  as  how  many  different  length  to  thickness  and  a/b 
ratios  may  be  investigated  for  the  various  y  directed  boundaries.  In  the 
subsequent  presentation  of  results  for  each  B.C. ,  the  range  of  S  and  plate 
aspect  ratio  for  which  useful  data  is  obtained  is  given.  (This  gives  a  good 
indication  of  the  limitations  of  the  program  written  for  this  thesis) .  For 
all  boundaries  however,  an  increase  in  the  length  to  thickness  ratios 
indicates  a  decrease  in  the  effect  of  shear  force  variations  through  the 
thickness.  That  is,  the  approximation  that  shear  reformation  effects  are 
negligible  is  increasingly  valid  as  the  plate  becomes  thinner  and  physical 
differences  in  the  z  direction  become  very  snail.  The  thinner  plates  quickly 


tend  to  the  classical  plate  behaviour  as  S  is  increased  past  40.  This  trend 

is  apparent  on  all  of  the  subsequent  figures  depicting  non-dimensional i zed 

2  2 

buckling  load  (-N^a  /*•  vs  length  to  thickness  ratio. 

The  non-dimens ionalizing  parameter  is  chosen  as  it  clearly  illustrates 
the  effects  of  two  important  variables  on  the  buckling  load.  The  first  is  the 
length  to  width  ratio  (S)  of  the  plate,  the  second,  the  aspect  ratio  of  the 
plate.  The  curves  assymptotically  approach  limiting  value  as  S  increases  end 
this  value  can  be  seen  to  represent  the  classical  laminated  solution.  No 
results  are  presented  for  length  to  thickness  ratios  under  10  for  any  of  the 
boundary  conditions.  It  is  important  to  realize  that  the  assumptions  of  plane 
stress  and  no  strain  in  the  z  direction  are  less  and  less  valid  as  S 


decreases.  The  errors  resulting  from  the  assumptions  are  no  longer  negligible 
at  length  to  thickness  ratio  less  than  10  and  any  results  obtained  are 
invalid.  By  choosing  this  non-dimensional izing  parameter,  the  reader  also 
gains  a  better  appreciation  of  the  dependence  of  the  buckling  load  on  the 
aspect  ratio  as  the  latter  is  varied  from  two  to  four. 


For  the  plate  simply-supported  on  all  sides,  a  square  configuration  is 
used  in  calculating  the  buckling  load.  Length  to  thickness  is  varied  frcm  io 
to  30  and  the  results  are  presented  in  Figure  3.1  for  the  first  bending  mode. 
The  curve  for  the  non-dimensionalized  buckling  load  is  seen  to  flatten  out  for 
S  approaching  32.5.  The  difference  at  this  point  between  CFT  and  shear 
deformation  theory  (SDT)  is  negligible.  For  a  thicker  plate  however,  the 
difference  increases  to  a  maximum  of  26.8%  for  an  S  of  10. 


For  the  plate  which  is  clamped  at  y=0  and  y=b,  a  rectangular  plate  with 
aspect  ratios  (a/b)  of  two,  three  and  four  is  used  in  the  analysis.  Length  to 


thickness  is  varied  from  10  to  37.5  for  the  first  plate.  The  curve  of  the 
non-dimensional ized  buckling  load  once  again  behaves  asymptotically  as  the 
value  of  S  becomes  greater  than  30  (See  Figure  3.2) .  For  an  aspect  ratio  of 
two,  a  maximum  difference  of  37.14%  exists  between  the  SOT  value  and  CPT 
extrapolated  value  at  S  equal  to  10.  As  the  aspect  ratio  increases,  the  non¬ 
dimens  ionl  ized  buckling  load  increases.  This  is  misleading  as  the 
dimensional  ized  buckling  load,  N^,  does  in  fact  decrease,  as  a/b  increases, 
indicating  a  decrease  in  the  plates'  ability  to  withstand  the  uniaxial 
compressive  force.  The  maximum  difference  between  CPT  and  SOT,  for  a/b  equal 
to  three,  is  approximately  36.65%.  For  a/b  of  four,  S  is  varied  from  20  to  40 
with  a  maximum  difference,  between  the  theories,  of  36.28%  occur ing  at  S  equal 
to  10.  CPT  solutions  for  orthotropic  plate,  presented  in  reference  [28], 
indicate  that  a  common  solution  exists  for  all  aspect  ratios  greater  than  two. 
Results  shewn  in  Figure  3.2  indicate  that  plates  of  different  aspect  ratios  do 
not  have  common  asymptotic  values  for  S  equal  to  40.  It  is  apparent  that  as 
a/b  is  increased  the  difference  between  asymptotic  values  does  not  vanish  due 
to  the  interaction  between  boundary  and  buckling  load. 

A  square  plate  is  used  to  analyze  the  simple-clamped  boundary  conditions 
in  y.  The  curve  in  Figure  3.3  is  plotted  for  length  to  thickness  ratios  of  10 
to  25.  Computer  limitations  prevent  larger  S  values  from  being  used  to 
calculate  the  budding  load.  The  maximum  difference  between  SOT  and  CFT 
occurs  at  a  length  to  thickness  ratio  of  10  and  has  a  magnitude  of  29.03%. 

The  extrapolated  value  for  the  classically  derived  buckling  load  is  reached 
fairly  quickly  at  S  equal  to  32.5. 


For  the  simple-free  boundaries  in  y,  the  behaviour  of  the  square  plate 
is  depicted  in  Figure  3.4.  S  is  varied  frcm  10  to  22.5,  the  upper  limit  a 
function  of  computer  limitations.  The  results  allow  the  characteristics  of 
this  boundary  to  be  studied  and  a  difference,  between  the  two  theories,  of 
23.01%  exists  at  S  equal  to  10  and  decreases  as  S  increases. 

The  clamped-free  y  directed  boundaries  are  studied  using  a  rectangular 
plate  with  aspect  ratios  varying  from  two  to  four.  For  a/b  of  two  and  three, 

S  is  increased  from  10  to  30.  For  a/b  of  four,  S  is  increased  to  40.  The 
three  curves  are  plotted  in  Figure  3.5.  Once  again,  though  the  curves  show 
asymptotic  behaviour  as  S  is  increased,  the  values  at  a  length  to  thickness 
ratio  of  40  are  different  for  each  plate  aspect  ratio.  This  demonstrates  the 
significance  of  the  shear  deformation  in  an  accurate  model  of  plate  behaviour. 
Maximum  differences  between  classical  theory  and  shear  deformation  theory 
varies  from  26.86%  to  27.84%  as  the  aspect  ratio  increases  from  two  to  four. 

The  final  results  obtained  for  a  buckling  problem  are  presented  in 
Figure  3.6.  The  free-free  boundaries  considered  are  analyzed  using  a 
rectangular  plate  with  three  different  aspect  ratios.  For  a  plate  aspect 
ratio  of  two,  the  length  to  thickness  ratio  is  varied  frcm  10  to  40.  A 
difference  of  26.92%  exists  between  SOT  and  CFT  at  S  equal  to  10  for  this 
first  rectangular  plate.  For  a/b  of  three,  S  is  also  varied  frcm  10  to  40  and 
the  disparity  between  theories  peaks  at  23.4%  for  S  equal  to  10.  The  final 
plate  studied,  with  an  aspect  ratio  of  four,  has  a  range  for  S  of  10  to  40. 

The  max  divergence  between  theories,  at  S  equal  to  10,  is  23.81%. 


s 


Figure  3.4  Plot  of  Normalized  Buck’ing  hoad  vs  Thickness 
Ratio  For  Simple-Free  Boundaries. 
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It  is  interesting,  at  this  point,  to  ocnpare  sane  of  the  results 

jjjj 

obtained  in  order  to  better  understand  the 

interaction,  if  any,  between  the 

M  jrii 

boundaries 

in  the  y  direction  and  shear  deformation  effects.  To  do  so,  plates 

with  similar  geometries  and  different  boundary  conditions  sire  compared  in 

1 

Table  3.5. 

i 

Boundaries  Aspect  Ratio 

Max  Difference 

K 

in  y 

of  Plate 

between  SDT/CPT 

1 

(%) 

| 

SS 

1 

26.80 

I 

SC 

1 

29.03 

SF 

1 

23.01 

CC 

2 

37.14 

| 

CC 

3 

36.65 

1  C# 

OC 

4 

28.00 

CF 

2 

26.86 

1 

CF 

3 

24.31 

; 

CF 

4 

27.84 

FF 

2 

26.92 

k> 

<N 

K 

FF 

3 

23.40 

ft 

FF 

4 

23.81 

1 

Table  3 

.5  Ccrrparison  of  Discrepancies  Between  Shear 

$ 

Deformation  Theory  and  Classical  Plate 

J 

f 

Theory  for  Stability  Problems. 

K 

g 

General  behaviour  characteristics  become  more  apparent  from  this  table. 

| 

First,  the 

effects  of  shear  deformation  are  very  important  regardless  of  the 

type  of  boundaries  for  the  plate.  An  average  difference  of  28.50%  exists, 
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regardless  of  the  B.C.  's  in  y,  between  the  models  of  plate  behaviour  for  the 
orientation  studied.  The  effects  of  shear  deformation  are  somewhat  lower  for 
the  free-free  boundaries  than  those  which  have  any  form  of  restraint.  The 
free-free  boundaries  have  an  average  difference  of  24.71%  compared  to  33.18% 
for  the  other  boundaries  in  y.  This  is  most  likely  due  to  a  greater  reaction 
between  bucJding  load  and  boundaries  if  the  boundaries  are  prevented  from  free 
movement. 

Secondly,  as  the  plates  get  thinner,  the  shear  deformation  effects 
remain  important.  For  the  rectangular  CF  or  FF  plates  under  study,  the  plate 
behaviour  becomes  less  influenced  by  the  shear  deformation  as  the  aspect  ratio 
is  increased.  In  addition,  the  asymptotic  value  reached  by  the  plate  as  it 
becomes  thinner  remains  different  for  the  range  of  aspect  ratios  investigated. 
Unlike  classical  theory  which  predicts  a  common  asyirptotic  value  for  all  a/b 
ratios  over  two  (refence  [28]) ,  SOT  indicates  that  shear  effects,  though 
small,  still  linger  even  as  S  increases  to  40.  Differences  in  the  asyirptotic 
value  do  decrease  however,  as  a/b  increases  for  the  free  boundaries. 

Thirdly,  the  results  validate  the  approach  presented  in  the  theory 
section  for  dealing  with  complex  conjugate  roots  to  Eq.  (65) .  These  roots 
occur  at  values  of  which  are  lower  than  for  the  case  of  three  real  roots. 
Hence,  an  erroneous  root  would  appear  higher  than  the  general  trend 
established  with  the  best-fit  curve.  For  the  stability  problem,  only  two  such 
possible  cases  did  occur  for  all  the  computations  and  they  are  indicated  on 
the  respective  figures. 


Fourthy,  the  plate  aspect  ratio  is  an  important  parameter  to  consider  in 
a  stability  problem  if  the  plate  boundaries  normal  to  the  applied  force 
provide  sane  restraint.  Figure  3.6  clearly  shows  that  free  boundaries  lead  to 
similar  values  for  critical  load  regardless  of  a/b. 

Finally,  stiffer  plate  supports  result  in  higher  buckling  loads.  This 
can  be  seen  by  comparing  values  of  N  calculated  for  the  CC  and  CF  boundaries. 


Vibration  Problem 


In  a  procedure  similar  to  the  one  followed  for  the  buckling  problem, 
each  different  boundary  condition  is  studied  by  elitering  the  length  to 
thickness  ratio  and  when  necessary,  the  length  to  width  ratio.  Initial  tests 
are  dene  with  rotary  inertia  effects  removed  from  the  model.  The  second  runs, 
with  RI  re-introduced  gives  a  relative  perspective  of  the  importance  of  shear 
end  inertia  in  an  accurate  representation  of  plate  action.  Generally 
speaking,  Table  3.4  provides  a  good  idea  of  the  trend  observed  for  sill  of  the 
vibration  problems  studied.  Rotary  inertia  is  found  to  decrease  the  overall 
stiffness  of  the  plate  by  an  average  value  of  less  than  one  percent.  Relative 
to  shear  deformation ,  rotary  inertia  complicates  a  plate  mathematical  model 
with  no  important  increase  in  accuracy. 

As  was  the  case  for  the  buckling  problem,  a  trend  that  is  observed, 

regardless  of  boundary  conditions,  is  the  diminished  inpact  of  shear 

deformation  as  plate  length  to  thickness  increases.  All  curves  of  u> 

2  2  1/2 

(ua  (p/E^h  )  )  vs  S  shew  a  decrease  of  the  discrepency  between  SOT  and  CPT 

as  S  increases.  Transverse  shear  variations  are  inverse  functions  of  plate 
thickness  and  become  negligible  for  high  enough  values  of  S. 

For  the  plate  sinply-supported  on  all  sides,  a  square  geometry  is  used 
in  calculating  the  natural  frequency.  Length  to  thickness  is  varied  from  15.0 
to  30.0  and  the  results  are  recorded  in  Figure  3.7.  The  curve  is  for  the 
model  that  included  rotary  inertia  effects.  The  average  difference 
with/without  RI  for  this  problem  is  0.19%.  A  difference  of  29.14%  exists 
between  SOT  and  the  extrapolated  CPT  when  the  length  to  thickness  is  10.  This 
difference  drops  very  quickly  as  the  curves  flattens  out  around  S  of  32.5. 


A  rectangular  plate  with  aspect  ratios  of  two,  three,  and  four  is  used 
in  the  analysis  of  the  clamped-clamped  y  boundary  conditions.  As  seen  in 
Figure  3.8,  S  is  varied  frcm  15  to  35  for  the  aspect  ratio  of  two.  Maximum 
deviation  between  theories  is  38.64%  for  this  geometry.  Variations  of  S  frcm 
20  to  35  are  used  in  calculations  for  a/b  equal  to  three.  It  becomes  clear 
from  Figure  3.8  that  a  narrower  plate  is  stiffer  as  the  natural  frequency 
increases.  The  plate  behaves  more  and  more  like  a  reinforced  beam  as  a/b 
increases  past  a  ratio  of  two.  Maximum  theoretical  difference  for  this 
configuration  is  33.4%  for  a/b  of  four,  this  maximum  readies  34.43%  as 
frequency  is  plotted  for  S  of  20  to  32.5. 

The  square  plate  is  once  again  used  to  analyze  the  simple-clamped 
boundary  in  y.  The  curve  in  Figure  3.9  is  valid  for  S  ranging  frcm  15  to  25. 
Computer  limitations  prevent  the  use  of  larger  S  values  in  the  computation  of 
the  natural  frequency.  A  maximum  divergence  of  36.4%  is  recorded  for  this 
plate  between  theory  considering  shear  through  the  thickness  and  classical 
analysis. 

For  the  simple-free  y  directed  boundaries,  plate  behaviour  is  depicted 
in  Figure  3.10.  A  square  and  a  rectangular  plate  are  used  to  show  frequency 
variations  as  functions  of  plate  length  to  thickness  ratios.  For  the  square 
plate,  with  a  variation  in  S  of  17.5  to  22.5,  the  maximum  difference  between 
theories  is  16.1%.  For  the  plate  with  aspect  ratio  of  two  and  variance  of  S 
of  15  to  35,  this  difference  peaks  at  20.2%.  For  this  problem,  it  is 
interesting  to  note  that  an  increase  of  a/b  frcm  one  to  two  represents  a 
decrease  in  plate  stiffness. 


The  clamped-free  plate,  whose  behaviour  is  illustrated  in  Figure  3.11, 
is  studied  for  three  different  aspect  ratios.  Once  again,  the  influence  of 
shear  deformation  is  very  noticeable  as  is  the  fact  that  this  influence  became 
negligible  quickly  as  S  approaches  40.  Maximum  differences  between  classical 
theory  and  shear  deformation  theory  vary  fran  19.44%  to  22.16%  for  aspect 
ratios  of  two  to  four  repectively . 

Important  trends  are  seen  to  be  slightly  different  for  the  vibration 
problem  than  those  recorded  for  the  buckling  problem.  First,  the  overall 
impact  of  shear  deformation  seems  to  be  equivalent  for  the  vibration  problem, 
as  illustrated  by  the  values  in  Table  3.6.  Regardless  of  y  boundary 
conditions ,  average,  difference  between  SOT  and  CPT  is  27.02%  for  the 
vibration  problem,  very  close  to  the  28.82%  obtained  for  the  buckling  problem. 
(In  this  comparison,  RI  effects  have  been  excluded  for  the  vibration  problem.) 
The  values  obtained  in  Table  3.6  compare  normalized  frequency  at  S  equal  to  10 
and  30.  For  values  of  S  less  than  10,  the  discrepancy  between  SOT  and  CPT 
quickly  increases. 

Secondly,  shear  is  more  important  to  consider  when  studying  the  boundary 
conditions  which  are  stiffer.  The  clairped-claitped,  simple-clamped  and  simple- 
simple  have  an  average  difference  between  the  two  theories  of  34.4%. 

The  boundaries  containing  a  free  edge,  on  the  other  hand,  have  an  average 
discrepancy  of  19.64%.  Thus,  the  shear  force  effects  are  higher  when  the 
influence  due  to  rigid  plate  boundaries  is  more  pronounced. 


Boundaries 
in  y 


Aspect  Ratio 
of  Plate 


Max  Difference 
between  SDT/CPT 
(%) 


29.14 


36.40 


16.10 


20.20 


38.64 


33.40 


34.43 


19.44 


20.32 


22.16 


Table  3.6  Oanparison  of  Discrepancies  Between  Shear 
Deformation  Theory  and  Classical  Plate 
Theory  for  Vibration  Problems 


Thirdly,  a  quick  ocrparison  shows  that  the  stiffness  of  the  boundaries 
does  affect  the  magnitude  of  the  natural  frequency  of  the  plate.  The  stiffer 
the  supports,  the  higher  the  natural  frequency  of  the  plate.  This  can  be 
clearly  seen  by  comparing  values  of  u>  for  the  CC  vs  the  CF  plate  or  the  SS  vs 
the  SF  plate. 


Fined  ly,  ocrparison  of  the  two  inch  plate  used  in  this  thesis  to  the  one 
inch  thick  plate  used  by  Bowlus  [2]  indicates  a  small  difference  in  the 
maximum  discrepancy  recorded  between  SDT  and  CPT.  This  is  expected  as  the 
thicker  heavier  plate  is  more  affected  by  force  variations  through  the 
thickness  as  they  tend  to  be  more  pronounced  than  for  a  thinner  plate. 


•'  r, 


Vl'oi' 


IV  Conclusions 


live  results  obtained  from  the  computations  performed  in  this  thesis 
allow  the  following  conclusions  to  be  presented.  They  include  comments  on  the 
Levy  technique  and  on  both  the  stability  and  vibration  problems.  It  should  be 
noted  that  all  conclusions  made  are  based  on  the  specific  laminate  used  in 
this  thesis.  Ihe  author  does  not  attempt  to  generalize  the  results  for 
laminates  of  arbitrary  compositions  nor  should  the  reader. 

Levy  Technique 

1.  Ihe  Levy  technique  is  a  viable  means  of  obtaining  base-line  solutions  for 
a  specific  class  of  laminated  ocnposite  plates. 

2.  Ihe  mathematics  of  the  solution,  especially  for  free  boundaries,  is  a 
great  deal  less  algabraically  complex  them  what  is  generated  using  a  Rayleigh- 
Ritz  or  Galerkin  approach. 

3.  Ihe  Levy  technique  cannot  be  extended  to  a  general  class  of  composite 
plates.  Ihe  presence  of  beniing-extensional ,  bending-twisting  coupling  terms 
or  D^/D.,  terms  would  not  allow  the  reduction  of  partial  differential 
equations  to  simple  differential  equations.  Thus,  though  the  Levy  method  has 
been  around  for  many  years,  it  has  never  been  fully  taken  advantage  of  due  to 
this  drawback. 

4.  Ihe  transcendental  equations,  which  are  solved  for  the  different 
boundaries  in  y,  are  very  sensitive  to  plate  geometry.  In  many  cases,  certain 
length  to  thickness  or  length  to  width  problems  could  not  be  resolved  due  to 
the  accuracy  limitations  of  the  computer  used. 


5.  Convergence  is  not  a  parameter  of  concern  with  the  Levy  procedure.  It  is 
not  a  term  dependent  approach  and  consequently,  the  accuracy  of  the  solution 
is  not  dependent  on  the  accuracy  of  the  displacement  models. 


Buckling  Problem 

1.  As  the  length  to  thickness  ratio  is  increased  past  40  for  the  plate,  the 
effects  of  shear  deformation  become  negligible.  Classical  plate  theory  can  be 
used  effectively  to  predict  plate  behaviour  for  the  thinner  plates. 

2.  Curves  of  the  nondimensionalized  buckling  load  vs  thickness  ratio  are  all 
mcmotonically  increasing  as  S  increases.  The  rate  of  increase  does  vary 
significantly  when  different  ranges  of  S  are  considered.  Increase  averages 
37.6%  for  S  from  10  to  20  and  28.7%  for  S  from  20  to  30  for  the  boundary 
conditions. 

3.  Shear  deformation  effects  account  for  an  average  difference  of  29.93% 
between  extrapolated  classical  theory  and  shear  theory.  The  effect  is  less 
pronounced  for  boundaries  conditions  which  are  not  very  stiff,  such  as  FF. 

Vibration  Problem 

1.  Curves  for  non-dimens iona  1  i zed  natural  frequency  vs  thickness  ratio  are 
all  monotonically  increasing  as  S  increases.  Rate  of  increase  does  vary  along 
the  curve,  being  the  largest  for  S  frcti  10  to  20  with  an  average  value  of 
38.8%. 

2.  Shear  deformation  effects  account  for  an  average  difference  of  27.02% 
between  extrapolated  classical  theory  and  shear  theory.  The  effect  is  more 


pronounced  for  the  boundaries  which  are  stiffer,  such  as  clanped-clarrped  or 
siirple-clanped . 

3.  Rotary  inertia  has  very  little  effect  on  the  overall  plate  stiffness  and 
can  be  neglected  in  the  mathematical  model  of  plate  behaviour  when  calculating 
the  first  natural  frequency.  For  higher  natural  frequencies,  RI  may  have 
significance  and  should  be  retained,  but  this  has  not  been  evaluated  within 
this  thesis. 
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Commands  used  to  obtain  simplified  expression  for  detCAijl 
1 J  determinant  ('/.)  ; 

2)  r  atexpand  ( '/.)  ; 

3)  xthru(X) ; 

4)  factor  (V.) ; 

5)  *  by  z**3/fhz; 


Expression  for  the  determinant  of  [Aij3. 
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4  2  2  2  2  2  2  2  2  2 

<  -  h  t  z  +  2fn  t  2  +  2  e  n  t  z  n  r2  -  c 

2  4  2  4  2  2  4  2  2  4 

-  g  h  n  t  2-fhn  t  2-f  n  t  2+e  n  t  z 

2  2  4  2  42  242 

-fn  rt  2+ghn  t  2  +  fgn  t  2  +  f  n  t  z 

4  £  2  4  4  2  : 

+  fn  r2-fgn  z  +  f  h  n  rt  -  g  h  n  rt  -  f 

£  £  4  2  4  2 

♦fgn  r+fht  2+gn  r2+e  n  rt)/(fhz) 


Expressions  for  Aj  in  Equat i on (£3) . 
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Expressions  for  B.j  in  Equation(81 ) . 

B3= f hz / f hz  =  1 . 0 
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Appendix  C 


Hie  first  program  is  a  fortran  rewrite  of  a  basic  program  used  by  Bowlus 
[2]  to  determine  the  nondimensional  stiffnesses  of  a  laminated  plate.  In 
particular,  the  program  is  rewritten  to  specifically  handle  [0/90]  ply 
layups  and  no  other  as  this  is  the  only  type  of  plate  studied  in  this  thesis. 
The  program  has  four  main  sections;  the  input,  the  ocnputation  of  lamina 
stiffnesses,  the  ocnputation  of  the  laminate  stiffnesses  and  the  output.  The 
input  obtains  the  following 

a)  plate  top  and  bottom  dimensions 

b)  E^,  E2,  G^,  and  mass  density  for  laminate  material 

c)  number  of  plies  in  laminate 

Based  on  this  information  and  the  fact  that  only  [0,90]^  laminates  are 
considered,  a  do-loop  is  enployed  to  calculate  the  non-dimensional  lamina 
stiffnesses  [Q^j  ]^.  An  "ms"  laminate  means  the  first  n/4  lamina  are  oriented 
at  0°,  the  next  n/2  at  90°  and  the  last  n/4  at  0°.  This  allows  the  use  of  a 
simplifying  do-loop.  Once  a  lamina  [Q. . is  calculated,  it  is  added  to  the 

lj  K 

sum  of  the  other  k-1  lamina  stiffnesses  and  the  k+1  stiffnesses  are 
calculated.  Doing  this  m  times  gives  the  non-dimensional  laminate  stiffnesses 
[Q^j  ] .  The  second  and  third  sections  of  the  program  completes  these  steps. 

The  output  is  presented  to  the  user  in  the  form  of  a  table  which  gives  the 
type  of  laminate  studied  and  a  listing  of  the  non-dimensional  stiffnesses 
calculated. 
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C «•*•**»****  ******************************************************** 

■:***  Program  to  compute  the  extensional  and  bending 
c *****  stiffness  elements  for  a  symmetric  laminate 
.;***♦*  build-up,  given  lamina  properties. 

c  *  ****************************************************************** 
wr i te<6, 10) 

10  for  mat  <’  enter  plate  top  i>.  bot  di  men  si  ons,  use  2fl0.5’,/> 

read (5, 15)  stop, z bot 
15  format <2f 10. 5) 

c 

c  section  to  obtain  lamina  data 
*2 

17  wr  i  t  e  <  6 ,  20 ) 

20  for mat (’ enter  El ,E2,G12, vl2  and  mass  density',/) 

wr i t  e ( 6 , 30 ) 

30  formatc'use  3el0.2  and  2f7.3,  OK?',/) 

r  ead  <  5,40  )  E 1 , E2 , 6 1 , V 1 , r ho 
40  format <3el0. 2, 2f 7. 3) 

45  write <6, 50) 

50  formatC'how  many  plies  in  this  plate, use  13',/) 

read (5, 60)  n 
60  format (13) 

V2=V1*(E2/E1> 

,G3=61 
62=0. 8*61 
al=0. 0 
a2=0. 0 
d  1=0.0 
d2»0.0 
d3=0. 0 
d4=0.0 
ps=0. 0 

C 

c  section  to  compute  A  and  D  elements  for  this  ply 

c 

Q 1 *E 1 / ( 1 . 0- C  V 1 * V2 ) ) 

Q2»  <  V 1 *E2 ) / < 1 . 0- <  V 1 *V2 ) ) 

Q3=E2/(1-(V1*V2) ) 

Q4-G2 

Q5»G3 

Q6-61 

t  k= (ztop-zbot ) /n 
k=n/4 

DO  100  1=1 ,n 

i f  < I . gt . k )  go  to  80 
70  th=0.0 

go  to  90 

80  if(I.gt.3*k>  go  to  70 

th«3. 1415927/2.0 
c 

c  section  to  compute  QBARS 

c  ! 

90  B 1=01 * (COS(TH) )**4+2. 0* (Q2+2. 0*Q6) *  < SINC  TH) ) **2* 

*  ( COS <  TH )  )**2-*-Q3*(SIN<TH)  )*#4 

B2= ( 0 1 +Q3-4 . 0*Q6 ) *  <  S I N  <  TH ) ) **2*  < COS ( TH ) ) **2+ 

*  Q2*C  CSIN(TH) )**4+<C0S<TH) )**4) 

B3=Q1 *  <  SIN  <  TH) ) **4*2. 0* (Q2+2. 0*Q6) * (SIN  <  TH) ) **2* 
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«. COS  <  TH  )  )  **2+Q3*  <: COS  <  TH  )  )  **4 

B4=  <  Q 1 +Q3-2 . 0*02-2 . 0*Q6 ) *  < S I N  < TH ) ) **2*  < COS  (  TH ) ) ** 
*  +Q6*< (SIN(TH) ) **4+ (COS(TH) )**4) 

B5=Q4#o:os<:th)  )**2+Q5*(SIN»:th)  :>**2 
B6=Q5*  < COS  <  TH )  )  **2+Q4*  (S  IN  <  TH)  >**2 

section  to  compute  lamina  A  elements 

Zl=ztop-I*tk 
ZK=z  t  op- ( I - 1 ) *t  k 
A3=B5* ( ZK-Z 1 ) 

A4=B6*( ZK-Z 1 ) 

F=RH0*«ZK-Z1> 


section  to  compute  lamina  D  elements 

DZ=<  ZK**3-Z 1**3) /3. 0 

F 1 =B 1 *DZ 

F2=B2*DZ 

F3=B3*DZ 

F4=B4*DZ 

section  to  find  laminate  A  and  D  elements 

A1=A1*A3 

A2=A2+A4 

D1=D1+F1 

D2-D2+F2 

D3«D3+F3 

D4-D4+F4 

PS»FS+P 

continue 

section  to  pri-nt  out  A  and  D  elements  of  laminate 
wr i te(6, 200) 

format  <26x,  ’graphite-epoxy  C  C> ,  90  ]  ’  ,  /,46x,  '2s’  ,  //, 
*32x,'one  inch  thick ’ , ///, 21 x, 9  element 4x , 

*’ dimensional  value’,///) 
wr i te<6, 300)  A1 , A2, Dl , D2, D3, D4 
format  <23x, ’ A44’ , 10x,F13.4, ///,23x, ’ A55’ , lOx, 
*F13.4,///,23x, ’Dll’ , lOx , F 1 3. 4, ///,23x, ’D12' , lOx, 

*F13. 4, ///, 23x, ’ D22’ , lOx , F 1 3. 4, / / / , 23x , ’ Dfcfc’ , 
*10x,F13.4,///, ’Units  for  Ai  .)  terms  are  lb/in’, 

*///, ’Units  for  Di j  terms  are  in-lbs’,//) 
wr i te (6, 350) 

for mat <  ’ another  problem  maybe?  yes= 1 , no®2’ , / ) 
read (5, 360)  j 
format (II) 

i f  <  j . eq . 2 )  go  to  400 
write(6, 370) 

format(’same  physical  pro.  but  diff.  #  plies?  yes=l’,/ 
read (5, 380)  m 
format (II) 
if(m.eq.l)  go  to  45 
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This  second  program,  also  written  in  Fortran,  has  four  main  sections. 

Die  transcendental  equations  in  the  fourth  section  of  this  program  are  solved 
using  the  incremental  search  technique.  A  flowchart  describing  the  mechanics 
of  this  technique  is  included  in  Appendix  E.  Die  four  sections  of  the  program 
are;  the  input,  the  calculation  of  the  coefficients  of  the  sixth  order 
equation  in  e,  Eq.  (63),  the  evaluation  of  a,  0,  7,  and  related  values  such  as 
A  (a) ,  and  the  solution  of  the  transcendental  equation  for  a  particular  set  of 
boundaries  in  y. 

Die  input  section  obtains  the  following  data: 

a)  number  of  modes,  n,  calculations  are  to  be  repeated  for 

b)  a/b  ratio  of  the  plate 

c)  does  user  want  to  include  effects  of  shear  deformation 
and/or  rotary  inertia  in  calculations? 

d)  the  boundary  conditions  in  y 

e)  buckling  or  vibration  problem  ? 

Die  following  are  an  intrinsic  part  of  the  program  but  can  of  course  be 
altered  by  editing  the  source  code: 

a)  the  material  properties  of  the  plate  under  investigation 

b)  a/h  ratio  used  in  calculation  (automatically  varied  frcm 
10  to  50  for  each  a/b  ratio  inputed) 

c)  initial  and  final  values  of  eigenvalue,  initial  increment 
for  eigenvalue  and  answer  tolerance  (all  parameters 
required  when  using  incremental  search) 


Depending  on  the  type  of  problem  being  solved,  the  second  section  of  the 
program  calculates  the  coefficients  of  Eq.  (63) .  Subroutines  are  used  for 
this  purpose  and  appropriate  terms  are  emitted  if,  in  the  vibration  analysis, 
rotary  inertia  is  removed.  Values  returned  are  used  in  the  third  section  to 
calculate  a,  fi,  and  7,  and  all  other  terras  required  to  define  the 
displacements.  These  terms  are  defined  by  Eqs.  (69),  (70),  (72),  (73),  (82), 
(83),  (84),  and  (85). 

The  expressions  for  the  transcendental  equations,  computed  in  the 
previous  section  of  the  thesis  and  subsequently  encoded,  are  evaluated  in  the 
fined  section  of  the  program.  An  eigenvalue,  for  the  particular  problem  under 
consideration,  is  obtained  when  the  transcendental  expression  is  identically 
zero.  In  using  this  program,  one  should  realize  that  certain  limitations  to 
the  plate  geometries,  which  may  be  studied,  do  exist.  These  limitations  are 
either  computational  or  theoretical  in  origin.  For  each  boundary  condition  in 
y  studied  in  the  thesis,  the  range  of  S  used  in  the  calculations  and  the 
aspect  ratio  of  the  plate  are  specified.  These  provide  the  user  with  a  very 
good  approximation  of  the  useful  oonputational  limits  of  this  program. 

Inconsistancies  present  in  the  theory  can  best  be  identified  by 
examining  the  graphical  results  and  observing  deviations  from  the  general 
trend.  As  was  discussed  earlier,  some  deviations  may  result  from  the 
assumption  that  all  roots  to  Eq.  (65)  are  real.  Neglecting  the  complex 
oonjugae  roots  is  a  good  simplification  as  only  two  instances  occur,  in  all 
computation,  where  the  roots  appear  to  be  complex.  A  second  source  of 
deviation  occurs  when  a  very  thick  plate,  S  around  ten,  is  studied.  In  this 
case,  the  assumptions  of  plane  stress  and  no  strain  in  the  z  direction  may 
lead  to  errors  which  are  no  longer  negligible,  depending  on  the  boundary 


conditions .  Such  deviations,  did  oocur  twice,  both  times  for  S  less  than  15 
and  are  identified  on  the  appropriate  figures.  Thus,  the  theory  is  inexact 
but  more  than  adequate  for  the  problems  studied  in  this  thesis. 
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c  #####*#*#####*#####***•*#***#*•*****#**•*******************#*#*#**#< 
c  Program  to  obtain  uniaxial  buckling  load  for  a  square 
c  composite  plate  given  boundary  conditions.  Program  can 
c  also  calculate  the  natural  frequencies  of  such  a  plate, 
c  ********************************* **##***#***##*#**#*##**#**##***#**^ 
DIMENSION  det  er ( 2 ) , A ( 3 ) , K APA (  3  > , KAP A 1(3), OMEGA (  3  >  ,  OMEGA 1(3) 

REAL  NU,  n 

INTEGER  type, be , pp, sd 

DOUBLE  PRECISION  aa, thi ck , r , u, coef t4, coef t2, coef tO, e, f , g , 

*h , z,x,b,c,d,p,q, ALPHA, BETA, GAMMA, A, KAPA, bb , vl , v2, v3, v4, 

*v5, v6, v7, vB, pi , BIGB, U, VE, ALPHA2, BETA2, GAMMA2, KAPA1 , 

#det , DIS, rho, al 1 , al2, al3, a21 , a22, a23, a31 , a32, a33, OMEGA,  1 

♦OMEGA 1 , eigen, xmax, dxi , epsi , del t x, v9, vlO, vll,vl2,vl3,vl4,o 
*E2, answer 

COMMON  /A/A44, D12, D66, Dll, D22, m, aa 
c  Material  properties  for  a  composite  plate  2"  thick. 

4  A44= 1080000. 625 

D1 1=12441315. 0 
D 12=281 690. 125 
D22=2582 159.0 
D66=400000. 0 
rho=0.055 
E2=l . 40d06 
'wr ite<6,6) 

6  format  (’  remember ,  thi  s  plate  is  2"  thick!’,/) 

10  write(6,40) 

40  format (’how  many  modes  problem  to  be  solved  for?, use  II’)  j 

read <5, 50)  mm 

c  The  user  must  choose  for  how  many  modes  he  wants  to  solve 
c  a  particular  a/h  &  a/b  geometry  for. 

50  format  (ID 

wr i te<6, 60) 

60  format (’problem  to  be  sol ved?buckl ing=l , vibs=2' ) 
read  (5,  70)'  type 
70  format  (ID 

wr ite(6,S0) 

80  format (’what  is  a/b  fc  thickness  of  the  plate?’) 
wr i te(6, 90) 

90  formatCute  2f8.3  to  input  val ues, okay?’ ) 

read(5,100)  ab, thick 
100  format (2f 8. 3) 
wr ite(6, 103) 

103  format(’do  you  want  to  include  SD?  yes=l,no*2’) 
read (5, 106)  sd 
106  format  (ID 

wr ite(6, 1 10) 

110  format(’do  you  want  to  include  RI?  yes=l,  no=2’ ) 
read(5,120)  L 
120  format  (ID 

wr ite(6, 125) 

125  format (’what  BC ’s  do  you  want  for  y=0,b  ?’,/) 
wr i te(6, 130) 

130  format (’SS- l, CC=2, SC-3, SF =4, CF=5,FF=6» , /) 

read (5, 135)  be 
135  format  (ID 

aaa*10.0 


c  Program  loops  through  m  from  1  to  mmy  the  number  specified 
c  earlier  by  the  user. 
aa-aaa*thick 
bb-aa/ab 

c  program  varies  a/h  from  10  to  55  for  every  value  of  a/b 
c  specified  by  the  user. 

145  k-1 

ei gen-200. 0 
xma  x -2000000 . 0 
dxi -100. 0 
epsi -0.00001 
deltx-dxi 

146  i f (type.eq. 1 )  go  to  150 
w-eigen 

go  to  160 
150  r -eigen 

c  ******************************************************* 
c  Start  of  second  section  of  program, 
c  Calculation  of  coefficients  of  angle  theta  used  to 
c  define  displacement  and  rotations. 

c  a###*#***#**#*#*#**# *********************************** 
call  buckle(coeft4, coeft2, coeftO, r ,n, e, f , g, h, z , sd) 
ifCsd.eq.2)  go  to  899 
'go  to  170 

160  call  vibs(coeft4,coeft2, coef to, wv rho, thick, n, e, f ,g, 
*h,L,x,z,sd) 
if(sd.eq.2)  go  to  899 

c  ****************************************************** 
c  Start  of  third  section  of  program. 

c  Evaluation  of  alpha,  beta,  gamma  and  related  terms, 

c  a***************************************************** 

170  b=coeft4 
c-coeft2 
d-coeftO 
p-c-b*b/3.0 

q-d-b*c/3.0+(b**3*2. 0/27.0) 

DIS— 27. 0*q*q-4. 0*p**3 
i f (DIS)300, 400, 400 

c  if  the  discriminant  is  <  0,  we  have  only  one  real  root 
c  to  the  cubic.  Physically,  this  does  not  make  sense.  The 
c  program  will  not  use  the  eigenvalue  in  any  further  calculations 
c  but  will  go  to  line  755.  Here,  the  eigenvalue  will  be  increased 
c  by  an  amount  of  deltx  and  b,c,and,d  will  be  recalculated. 

300  go  to  755 
400  pi -3. 1415927 

do  500  1-1,3 

BIGB-dsqrt (-4. *p/3. ) 
ll-3.*sqrt  (3.  )*q/<2.  *p*dsqrt  C-p)  ) 

VE= (dacos (U) /3. )  +  < 1-1 >*pi*2.  /3. 

NU-dcos(VE) 

A(I)=NU*BIQB-b/3. 

500  continue 

ALPHA2=A<1) 

BETA2-A<2)  # 

GAMMA2=A(3) 

i  f  (  ALPHA2 . 1 1 . 0 .  )  ALPHA2—  1 .  * ALPHA2 
i  f  (BETA2.  lt.0.  )  BETA2—  1 .  *BETA2 
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i f  <  GAMMA2 .  1 1 . 0 .  >  GAMMA2=- 1 . *GAMMA2 
ALPHA=*dsqrt  (ALPHA2) 

BETA*dsqrt (BETA2) 

GAMMA=dsqrt (6AMMA2) 
i f (type. eq. 2)  go  to  650 
do  550  J«l,3 

c  Buckling  Problem 

c  if  one  of  the  roots  of  the  cubic  is  negative,  then  eqsdJU 

c  and  (73)  must  be  used  instead  of  eqs(69)  and  (.TO). 

if (A(J).lt.O.  >  go  to  530 
c  eq(fc9)  and  eq(lo)  follow: 

KAPA< J>=( (h*A(J)-z-n*n#f )*n*z-n*e*z*A( J) )  . 

*  /  ( <h*A<  J)-z-n*n*f  )*<  f*A(  J)-z-n*n*g)+ 

*  n*n*e*e*A( J) ) 
if(J.eq.l)  angle= ALPHA 
if(J.eq.2)  angle=BETA 
if(J.eq.3)  angle-GAMMA 

KAPA1  ( J)=(n*z-(  f*A(  J)-z-n*n*g)*KAPA(  J) )  /  (n*e*angle)  £fiiJ.(1o) 
go  to  550 

c  eq(7>,  and  eq(73)  follow: 

530  KAPA( J>=( (h*A( J)+n*n*f+z )*n*z-n*e*z*A( J) ) 

*  / (n*n*e*e*A( J)-(h*A< J)+n*n*f+z )*( f*A( J)+n*n*g+z ) )  W*'“  ' 

ifCJ.eq.l)  angle=ALPHA 

'if(J.eq.2>  angle=BETA 
if(J.eq.3)  angle=GAMMA 

KAPA1 C  J)=(n*z+(  f*A<  J)+n*n*g+z  )*KAPA(  J) )  /  (n*e*angl e)  £fe|0>tl3) 
550  continue 

go  to  700 

650  do  685  JJ-1,3 

c  Vibration  Problem 

c  if  one  of  the  roots  of  the  cubic  is  negative,  then  eqs(P*f) 

c  and  (  0*  >  must  be  used  instead  of  eqs(&2)  and  (©3). 

i f (A( JJ) . 1 t . 0. )  go  to  660 
c  eq(fX)  and  eq(f3)  follow: 

0MEGA( JJ)=( (h*A( JJ)-z-n*n*f+x)*n*z-n*e*z*A(JJ) ) 

*  /( (h*A(JJ)-z-n*n*f+x)*<  f*A(JJ)-z-n*n*g+x)+  r 

*  n*n*e*e*A( JJ) > 
if(JJ.eq.l)  angle=ALPHA 
ifCJJ.eq.2)  angle=BETA 
if(JJ.eq.3)  angle=GAMMA 

OMEGA  1  (JJ)=(n*z-(  f*A( JJ)-z-n*n*g+x)*0MEGA( JJ)  )/ 

*  (n*e*angle) 
go  to  685 

c  eq ( M )  and  eq(l$)  follow: 

660  OMEGA ( JJ>=( (h*A( JJ)+n*n*f+z-x)*n*z-n*e*z*A( JJ) ) 

*  / (n*n*e*e*A( JJ)-(h*A( JJ)+n*n*f+z-x)*( f*A( JJ)+n*n*g+z-x) > 

if(JJ.eq.l)  angle=ALPHA  mMBH) 

if(JJ.eq.2)  angle=BETA 

ifCJJ.eq.3)  angle=GAMMA 

OMEGA  1  (JJ)=(n*z+(  f *A ( JJ) +n*n*g+z-x ) *0MEGA ( J J) ) /  (gs) 

*  (n*e*angle) 

685  continue 

c  ****************************************************************** 
c  Start  of  fourth  section  of  program. 

c  solution  to  transcendental  equation  of  boundary-value  problem, 

c  ****** ************************************************************ 

700  goto  (960,1000,1100,1200,1300,1400)  be 
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•  deter (k)*det 

if(k.eq.2)  go  to  790 
k-k+1 

>  eigen®eigen+deltx 
go  to  146 

'  i f (deter (l)*deter (2) )  800,820,795 
i f (eigen. gt . xmax)  go  to  850 
deter ( 1 )=deter (2) 
go  to  755 

>  if (deltx-epsi)  820,820,810 
ei gen=ei gen-del t x 

delt x=del tx/10. 
go  to  755 

i f (type. eq. 2)  go  to  825 
answer=eigen*bb*bb/(pi*pi*aa*sqrt (g*h) ) 
go  to  827 
answer =ei gen 

write(6,830)  aa, bb , m, answer 

format ( ' aa®’ , f8. 3, x, 'bb®' , f8. 3, x, 'n*' , II, x, ' root®' , d20. 10, / 
go  to  890 
wr ite(6,860) 

format ('max  value  reached  and  no  root  found',/) 
if(m.ge.mm)  go  to  895 
'm=m+ 1 
go  to  145 

i f (aaa. ge. 45. 0)  go  to  899 
aaa®aaa+2.5 
go  to  140 
wr i t  e ( 6 , 900 ) 

format ('do  you  want  to  try  another  probl em?yes=l ' ) 
read (5, 910) J 
format  (ID 

i f ( j.ne. 1)  go  to  950 
wr i te(6, 920) 

format(*any  change  in  mat’l  properties?  y=l',/> 

read (5, 930)  pp 

format  (ID 

if(pp.eq.l)  go  to  4 

go  to  10 

stop 

i f (a (2) . gt . 0. 0)  go  to  970 

Routine  to  calculate  detCaijl  for  "simple-simple"  be  in  y. 
det=beta-pi /bb 
go  to  750 

i f (type. eq . 2)  go  to  980 

vl®kapa( 1) 

v2®kapa(2) 

v3®kapa(3) 

v4*kapal (1) 

v5»kapal (2) 

v£=kapal (3) 

go  to  990 

vl®omega( 1 ) 

v2«omega(2) 

v3®omega(3) 

v4®omegal ( 1 ) 

v5®omegal (2) 


v6somegal (3) 

990  al 1= Cdsi nh (al pha*bb) ) / 1 . dlO 

al2=(dsinh(beta*bb) )/l.dlO 
al3=(dsinh (gamma*bb) ) /l.dlO 
a21=(vl*dsinh (alpha*bb) )/ l.dlO 
a22=(v2*dsinh(beta*bb) )/l.dlO 
a23= ( v3*dsi nh (gamma*bb ) ) / 1 . d 10 
a31= ( v4*al pha*dsi nh (al pha*bb) ) / 1 . dlO 
a32= ( v5*bet a*dsi nh (bet a*bb ) >/l .dlO 
a33=( v6*gamma*dsi nh (gamma*bb) ) /I . dlO 
det=al I*(a22*a33-a32*a23>-al2*(a21*a33-a3l*a23)+ 

*al3* (a21*a32-a31*a22) 
go  to  750 
stop 

1000  i f (type. eq. 2)  go  to  1010 

c  Routine  to  calculate  detCAij]  for  "clamped-c lamped"  be  in  y. 
vl=KAPA(l> 
v2=KAPA(2) 
v3=KAPA(3) 
v4=KAPAl ( 1 ) 
v5=KAPAl (2> 
v6=KAPA 1 (3) 
v7= ( vl-v3) / ( v3-v2) 

'  v8= ( v2- v 1 ) / ( v3-v2 ) 
go  to  1020 

1010  v 1 =0ME6A ( 1 ) 
v2=0ME6A(2) 
v3=0MEGA (3) 
v4=0MEGAl (1) 
v5=0MEGA 1(2) 
v6=0MEGAl (3) 
v7=(vl-v3)/(v3-v2> 
v8=*  (  v2- v  1 )  /  <  v3- v2  ) 

1 020  a 1 1 *  <  ds i nh  <  ALPHA*bb ) - v4/ v6*ds i nh ( GAMMA*bb >  >  / 1 .  d  1 0 
a 1 2= ( dc osh  <  ALPHA*b  b ) + v7*dc  os ( BET A*bb ) 

*+ v8*dc osh (GAMMA*bb) >/ l.dlO 
al3=(dsin (BETA*bb)+v5/v6*dsinh (GAMMA*bb) ) /I . dlO 
a2 1 * ( v 1 *dsi nh (ALPHA*bb ) - v3*v4/ v6*dsi nh (GAMMA*bb ) ) / 1 . d 1 0 
a22=(vl*dcosh(ALPHA#bb>+v2#v7*dcos(BETA*bb> 

*+ v3* v8*dc  osh (GAMMA*bb ) ) / 1 . d 10 
a23= ( v2*dsi n ( BETA*bb ) +v3*v5/v6*dsi nh (GAMMA*bb > ) / 1 . d 10 
a31  =  ( v4* ( dc  osh ( ALPH A*bb ) -d  c  osh ( GAMMA*bb ) ) ) / 1 . d 1 0 
a32= ( v4*dsi nh ( ALPHA*bb ) +v5*v7*dsi n ( BETA*bb ) 

*+vS*v8*dsinh (GAMMA*bb ) ) / 1 . d 10 
a33= <  v5* ( -dc  os ( BETA*bb ) +dc  osh ( GAMMA*bb ) ) ) / 1 . d 1 0 
det~al I*(a22*a33-a32*a23)-ai2*(a21*a33-a31*a23>+ 
*al3*(a21*a32-a31*a22) 
go  to  750 
stop 

1100  i f (type.eq. 2)  go  to  1110 

c  Routine  to  calculate  detCAij]  for  "simple-clamped'  bc's  in  y. 
v 1 =KAPA ( 1 ) 
v2=kapa(2> 
v3=kapa(3) 
v4»kapal ( 1 ) 
v5=kapal (2> 
v6»kapal (3) 


go  to  1120 

1110  vl=omega(l) 
v2=omega(2) 
v3=omega(3> 
v4=omegal ( 1 ) 
v5=omegal (2) 
vG=omegal (3) 

1120  if <a(2).lt.0.0>  go  to  1130 
v20=dsi nh  <bet  a*bb > 
v30=dcosh (beta*bb ) 
go  to  1140 

1130  v20=dsin (beta*bb) 
v30=-dc  os  <  bet  a*bb ) 

1140  al 1= (dsi nh (al pha*bb ) ) / 1 . dlO 
a 1 2s  <  v20 ) / 1 . d 1 0 
a 13= (dsi nh (gamma»bb) ) /I .dlO 
a21=(vl*dsinh(al pha*bb ) ) /l . dlO 
a22=(v2*v20) / 1 . d 10 
a23= ( v3*dsi nh (gamma#bb ) ) / 1 . dlO 
a31=( v4*dcosh (alpha*bb) ) / 1 . dlO 
a32= <  v5* v30  > / 1 . d 1 0 
a33=C v6*dcosh (gamma*bb) ) / 1 . dlO 

det  =a 1 1 * ( a22*a33-a32*a23 ) -a 1 2*  <  a2 1 *a33-a3 1 *a23  >  + 

•a 1 3* ( a2 1 *a32-a3 1 *a22 ) 
go  to  750 
stop 

1200  i f (type.eq. 2)  go  to  1210 

c  Routine  to  calculate  detCAijl  for  "simple-free"  be  * s  in  y. 
vl=kapa( 1 ) 
v2=kapa(2) 
v3=kapaC3) 
v4=kapal < 1 > 
v5=kapal (2) 
v6=kapal (3) 
go  to  1220 

121C  vl=omega(l> 
v2=omega(2> 
v3=omega (3) 
v4=omegal C 1 ) 
vS=omegal <2) 
v6=omegal (3) 

1220  al 1=( (h*alpha*v4-n*vl*(e-f > ) *dsinh (alpha«bb> ) / 1 . dlO 
al2=( <h*beta*v5-n*v2*<e-f ) )*dsin<beta*bb) ) / 1 . dlO 
al3=( <h*gamma*v6-n*v3*  <e-f ) )*dsinh(gamma*bb) ) /I . dlO 
a21  =  ( <alpha*vl+n*v4)*dcosh  Caipha*bb> ) /I . dlO 
a22=( <beta*v2-n*v5)*dcos<beta*bb) )/l . dlO 
a23=( < gamma* v3+n«v6) *dcosh (gamma*bb) )/l.dlO 
a31=( < alpha* v4)*dcosh Calpha*bb) )/l .dlO 
a32=( <beta-v5) *dcos (bet a*bb ) ) / 1 . d 10 
a33= ( ( gamma* v6 ) *dc osh  <  gamma*bb ))/l.dlO 
det -a 1 1 • ( a22*a33-a32*a23 ) -a 1 2* ( a2 1 *a33-a3 1 *a23 >  * 
*al3«(a21*a32-a31*a22> 
go  to  750 
stop 

1300  i f (type. eq. 2)  go  to  1310 

c  Routine  to  calculate  detCAiJl  for  "c 1 amped-free"  bcrs  in  y. 
vl-kapa(l) 
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v2*kapa<2> 
v3>kapa(3) 
v4skapal ( 1 ) 
v5»kapal(2> 
v6-kapal(3> 
v7«(vl-v3)/(v3-v2> 
v8“< v2-vl ) / < v3-v2) 
go  t  o  1 320 

1310  vlsomega(l) 
v2=omega(2) 
v3=omega<3) 
v4=omegal ( 1 ) 
v5=omegal (2) 
v&=omegal (3) 
v7«(vl-v3)/(v3-v2> 
v8=(v2-vl ) / ( v3-v2> 

1320  al  ls( (-n*(e-f )*vl+h*alpha*v4)*dsinh(alpha*bb>+ 

*< v3*v4/v6*n*  <e-  f  )  -h*v4*gamma)*dsinh(gamma*bb)  )  / 1 .  dlO 
a 1 2*  <  C -n*  < e- f ) * v 1 +h *a 1 pha* v4 ) #dc osh  C al pha*bb )  + 

*(-n*(e-f )*v2+h*beta*v5)*v7*dcos(beta*bb)+ 

*  < -n*  <e- f ) *v3+h*gamma*v6 ) * v8*dcosh ( gamma*bb ) ) / 1 . d 10 
al3=( (-n*(e-f )*v2+h*beta*v5)*dsin (beta*bb)+ 

*(-n*(e-f )*v3*v5/ v6+h*gamma*v5> *dsinh( gamma *bb) ) /I .d 10 
'  a21  =  (  (alpha*vl+n*v4)*dcosh(alpha*bb>- 
*< gamma* v3* v4 / v6+ v4 ) *dc osh (gamma*bb ) ) / 1 . d 10 
a22=( (alpha*vl+n*v4)*dsinh(alpha*bb>-(beta*v2+n*v5)*v7* 
*ds i n (bet  a*bb )  ♦  < gamma* v3+n*v6 >  * v8*dsi nh ( gamma*bb >>/l.dl0 
a23= <  <  bet  a* v2-n* v5 ) *dc  os ( bet  a*bb  >  + 

*  (gamma*v3*v5/ v6+n*v5> *dc osh (gamma*bb ) ) / 1 . d 10 

a3 1 * ( ( al ph  a+ v4 ) *dc  osh ( a 1 pha*bb  > - <  gamma* v4/ v6+ v4 ) * 

*dcosh (gamma*bb ) > / 1 . d 10 
a32= ( (al pha+ v4 ) *dsi nh ( al pha*bb )  + (-bet  a+v5>  *v7* 

*dsin  (beta*bb)  +  (gamma-*-v6)*v8*dsinh(gamma*bb)  )  /I .  dlO 
a33=  (  (bet a'- v5 )  *dc os  ( bet a*bb  >  +  ( gamma* v5/v6+v5 )  • 

*dc  osh ( gamma*b  b > ) / 1 . d 1 0 

det  «a 1 1 * ( a22*a33-a32*a23 ) -a 1 2* ( a2 1 *a33-a3 1 *a23 >  + 

*a 1 3* ( a2 1 *a32-a3 1 *a22 > 
go  to  7S0 
stop 

1400  i f (type.eq. 2)  go  to  1410 

c  Routine  to  calculate  detCaij]  for  "free-free"  bc's  in  y. 
vl=kapa( 1 ) 
v2=kapa(2) 
v3=kapa<3) 
v4=kapal ( 1 ) 
v5-kapal (2) 
v6=kapal  (3) 
go  to  1420 

1410  vl*omega(l> 
v2=omega(2) 
v3*omega(3) 
v4somegal ( 1 ) 
v5»omegal (2) 
v6«omegal (3) 

1420  v7«-(n*v4+vl >/(v3+n*v6> 

if  (a(2>.gt.O.)  go  to  1430 
v8«(n*v5-v2>/ (n*v&+v3> 
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v9-  C  v4+al  pha )  +  <  v6+gamma  )  *  v7 
vl0*(beta-v5)+(v6+gamma)*v8 
go  to  1440 

1430  v8- (n*v5+v2) / (n*v6+v3) 

v9-  C  v4*al pha) ♦ <  vS+gamma ) *v7 
v 10* (beta*v5) ♦ < v6+gamma) *v8 

1440  vll—  v9/vl0 

vl2-v7+v8*vl 1 

vl3* (h*alpha*v4-n* (e-f )*vl ) / Cn* (e-f ) *v3-h*gamma*v6) 
vl4* (h*beta*v5-n* (e-f ) *v2) / (n* (e- f ) *v3-h*gamma*v6> 
o=e-f 

i f <a<2>. It. 0.0)  go  to  1450 

al l-( (-n*o*vl+h*v4*alpha)*dsinh(alpha*bb)+(-n*o*v2*vl 1 
*+h*v5*bet a*v 1 1 ) *dsi nh (bet a*bb ) ♦ <-n*o*v3*v 1 2+h* v6*v 1 2* 
•gamma) *dsinh (gamma*bb ) ) / 1 . d 10 
a 1 2* ( < -n*o*v 1 +h* v4*al pha ) *dc  osh ( al pha*bb )  + 
*(-n*o*v3*vl3+h*v6*vl3*gamma)*dcosh (gamma*bb) ) / 1 . dlO 
a 1 3=  <  < -n*o*v2+h* v5*bet a ) *dc osh (bet a«bb )  + 

* (-n*o*v3*vl4+h*v6*vl4*gamma) *dcosh (gamma*bb) ) / 1 . d 10 
a21* ( (vl*alpha+n*v4)*dcosh (al pha«bb )+ 

*(v2*beta*vl l+n*v3#vl l ) *dcosh (beta*bb ) + 

* ( v3*vl2*gamma+n*v6*vl2)*dcosh (gamma*bb) ) /I . dlO 
a22= ( (vl*alpha+n*v4) *dsinh (alpha*bb) ♦ 
*(v3*vl3*gamma+n*v6*vl3)*dsinh(gamma*bb) )/l .dlO 
a23= ( ( v2*beta+n*v5) *dsi nh (bet  a*bb ) ♦ 

* ( v3*vl4*gamma+n*vfc*vl4)*dsinh (gamma*bb) ) /I . dlO 
a31«( (v4+alpha)*dcosh (alpha*bb)+(v3*vl l+beta*vl 1 )* 

•dcosh (beta*bb )+( v6*vl2+gamma*vl2)*dcosh(gamma*bb) ) / 1 .dlO 
a32»(  (v4+alpha)*dslnh(alpha*bb)-«'(v6*vl3+gamma*vl3)* 
*dsinh(gamma*bb) )/l.dlO 

a33*«  (  ( v5+bet  a )  *ds i  nh  ( bet a*bb )  ♦  (  v6*v  1 4 +gamma* v  14)* 
•dsinh(gamma«bb) ) / 1 . dlO 
go  to  1460 

1450  al 1*( (-n*o*vl+h*v4*alpha)*dsinh(alpha*bb)+(-n*o*v2*vl 1 
*+h*v3*beta*vl 1 )*dsin(beta*bb)*(-n*o*v3*vl2+h*v6*vl2* 
•gamma) *dsi nh (gamma«bb ) ) /l . dlO 
a  12* ( (-n*o*vl+h*v4*alpha) *dcosh (alpha*bb)  + 
*(-n*o*v3*vl3+h*v6*vl3*gamma)*dcosh (gamma*bb) ) /I . dlO 
a  1 3* ( ( -n*o*v2+h*v5*bet  a ) *dc  os ( bet  a*bb ) ♦ 

* (-n*o*v3*vl4+h*v6*vl4*gamma) *dcosh (gamma«bb) ) /I . dlO 
a21*( ( vl*alpha+n*v4) *dcosh (alpha*bb)+ 

*(v2*beta*vl l-n*v5*vl 1 ) *dcos (beta*bb)+ 

*  ( v3*v  1 2*gamma-*-n*v6*v  1 2 )  *dc osh  (gamma *bb ) )  / 1 .  d  10 
a22* ( ( v 1 *al pha+n*v4 ) *dsi nh ( al pha*bb ) ♦ 

* ( v3*v 1 3*gamma+n* v6*v 13) *dsi nh ( gamma«bb ) ) / 1 . d 1 0 
a23=( (-v2*beta+n*v5)*dsi n(beta*bb)+ 

* ( v3*v 14*gamma+n*v6*vl4 ) *dsi nh (gamma*bb ))/l.dl0 
a31*( (v4+al pha) *dcosh (al pha*bb)* (-v5*vl l+beta*vl 1 )* 

•dc  os ( bet  a*bb )  +  ( v6*v 1 2+gamma*v 12) *dc  osh ( gamma*bb ) ) / 1 . d 10 
a32* ( (v4+alpha) *dsi nh (al pha*bb)+ (v6*vl3+gamma*vl3)* 
•dsinh (gamma*bb ) ) / 1 . d  10 

a33* ( ( v5-bet  a ) *dsi n (bet  a*bb ) ♦ ( v6*v 1 4+gamma*v 14)* 

•dsi nh (gamma *bb ) ) / 1 . d  10 

1 460  det -all * (a22*a33-a23*a32 ) -a 1 2* ( a2 1 *a33-a3 1 *a23 ) + 

•a 1 3* ( a2 1 *a32-a3 1 *a22 ) 
go  to  750 
stop 
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•nd 

SUBROUTINE  buckle(coef t4, coef t2, coef tO,  r f n, 
***ff9»h,z,sd> 
c 

c  Subroutine  to  calculate  coefficients  of  angle 
c  theta  ,  for  the  buckling  case 
c 

COMMON  /A/A44, D12t  D66, Dll, D22, m, aa 
REAL  n 
INTEGER  sd 

DOUBLE  PRECISION  coef t4, coef t2, coef tO, r , e, 

*f ,g,h,z, pi ,aa 
pi*3. 1415927 
if(sd.eq.l)  go  to  1910 
wr i te (6, 1900) 

1900  format< ’you  have  removed  shear  deformation  effects’) 
wr i te (6, 1901) 

1901  format ('from  the  problem.  You  no  longer  have  a  cubic’) 
wr  ite(6, 1902) 

1902  f ormat ( ’ equat i on  to  solve  so  that  the  rest  of  this’) 
wr  i te(6, 1903) 

1903  format (’ program  is  useless.  An  altogether  different’) 
wr i te(6, 1904) 

1904  ' format (’ program  would  have  to  be  used.  Sorry!!’,/) 

return 

1910  z=S./6.*A44 

e=D12+D66 
f«D66 
g-Dl  1 
h-D22 

n»FL0AT (m)*pi /aa 

coeft4=-l . *<h*z*z+g*h*n*n#z+f*h*n*n*z+f*f*n*n*z- 
*e*e*n*n*z-f *h*n*n*r ) / <  f *h*z ) 
c oef  1 2=*- 1 .  *  (  -2.  *  f  *n*n*z  *z -2.  *e*n*n*z  *z  +h*n*n*r  *z  ♦ 
#f*n*n*r *z-g*h*z*n**4-f*g*z*n**4-f*f*z*n**4+ 
*e*e*z*n**4+g*h*r *n**4+f*f *n*n*r-e*e*r*n**4) / <  f*h*z ) 
coeft0“-l . *<-n*n*z*z*r-g*z*r*n**4-f*z»r *n**4-f*g*r*n**6 
*+g*z*z*n**4+f *g*r *n**6) / <  f*h*z ) 
return 
end 

SUBROUTINE  vi bs (coef t4, coef t 2, coef tO, w, rho, thick, n, 

*e , f , g , h , L , x , z , sd ) 
c 

c  Subroutine  to  calculate  coefficients  of  angle 
c  theta,  for  the  vibration  case, 
c 

COMMON  /A/A44, D1 2, D66, D1 1 , D22, m, aa 
REAL  n 
INTEGER  sd 

DOUBLE  PRECISION  coef t 4, coef t 2, coef tO, w, rho, thick , 

•e, f,g,h,z,x,pi,aa 
i f ( sd . eq . 1 )  go  to  1 960 
wr i te<6, 1950) 

1950  format(’you  have  just  decided  to  remove  all  shear  thru’) 
write(6, 1951) 

1951  format (’the  thickness  effects  from  this  problem.  The') 
wr  i te(6, 1952) 
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1952  format ('program  you  arm  using  is  no  longmr  valid  as  it  is') 
wr i t  m ( 6 , 1953) 

1953  formatCsmt  up  to  solvm  a  cubic.  A  totally  diffmrmnt') 
wr i te(6, 195*) 

1954  format (' program  would  havm  to  bm  usmd  for  this  simplmr') 
wr i te(6, 1955) 

1955  format ('cast.  So  sorry  old  chap.',/) 
rmturn 

1960  z«5./6.*A44 

pi«3. 1415927 
e=D12+D6£ 
f«D66 
g*Dll 
h*D22 

if(L.aq.l)  go  to  2000 
x=0. 0 

go  to  2010 

2000  x=w*w*rho*thi ck**3/ ( 144. *32. 174) 

2010  y=rho*thick*w*w/ <12. *32. 174) 
n=FL0AT (m) *pi /aa 

cocf t4*(-h*z*z+h*x#z+f *x*z-g*h*n*n*z-f *h*n*n*z 
*-f*f*n*n*z+e*e*n*n*z+f*h*y)/< f*h*z ) 
coef t2= <-x*z*z+2.  *f*n*n*z*z+2.  *e*n*n*z*z-h*y#z-f*y*z 
*+x*x*z-h*n*n*x*z-g*n*n*x*z-2. *f *n*n*x*z+g*h*z*n**4 
*+  f *g*z  *n*#4+ f #  f *z*n**4-e*e*z  *n**4+h*x*y+  f *x*y 
*-g*h*n*n*y-f * f*n*n*y+e*e*n*n*y ) / ( 1 *h#z > 
comf t0*<y*z*z+n*n*x*z*z-g*z*z*n**4-2. *x*y*z+g*y*z*n*n 
*+f*n*n*y*z-n*n*x*x*z+g*x*z*n**4+f*x*z*n**4-f*g*z*n**6 
*+x*x*y-g*n*n*x*y-f *n*n*x*y+f *g*y*n**4) / ( f*h*z ) 
return 
end 
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19.  ABSTRACT 

'»  ® 

An  analytical  study  is  conducted  to  determine  the  stability  and  free 
vibration  characteristics  of  laminated  anisotropic  plates  using  the  Levy 
approach.  Included  in  the  plate  model  are  the  effects  of  shear  deformation 
and  rotary  inertia.  Six  different  boundary  conditions  in  the  y  direction  are 
analyzed  in  conjunction  with  simply-supported  boundaries  in  the  x  direction. 
The  y  directed  boundaries  considered  are  simple-simple,  clamped-clamped , 
simple-clamped,  simple-free,  clamped-f ree ,  and  free-free. 

Solutions  are  presented  for  the  buckling  loads  and  natural  frequencies 
of  rectangular,  graphite-epoxy  symmetric  plates.  The  results  indicate  the 
imprtance  of  including  shear  effects  and  rotary  inertia  in  a  plate's 
mathematical  model.  The  overall  importance  of  these  equation  parameters  is 
definitely  a  function  of  the  boundary  condition  and  a  general  statement  cannot 
be  made.  In  addition,  the  effectiveness  of  the  Levy  technique,  in  studying 
laminated  problems,  becomes  apparent  in  handling  the  more  complicated 
boundaries  as  compared  to  the  Galerkin  or  Rayleigh-Ritz  techniques. 


